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BREAKING  SPACES  AND  FORMS  FOR  THE  DPG  METHOD  AND 
APPLICATIONS  INCLUDING  MAXWELL  EQUATIONS 


C.  CARSTENSEN,  L.  DEMKOWICZ,  AND  J.  GOPALAKRISHNAN 

Abstract.  Discontinuous  Petrov  Galerkin  (DPG)  methods  are  made  easily  implementable 
using  “broken”  test  spaces,  i.e.,  spaces  of  functions  with  no  continuity  constraints  across 
mesh  element  interfaces.  Broken  spaces  derivable  from  a  standard  exact  sequence  of  first 
order  (unbroken)  Sobolev  spaces  are  of  particular  interest.  A  characterization  of  interface 
spaces  that  connect  the  broken  spaces  to  their  unbroken  counterparts  is  provided.  Stability 
of  certain  formulations  using  the  broken  spaces  can  be  derived  from  the  stability  of  analogues 
that  use  unbroken  spaces.  This  technique  is  used  to  provide  a  complete  error  analysis  of  DPG 
methods  for  Maxwell  equations  with  perfect  electric  boundary  conditions.  The  technique 
also  permits  considerable  simplifications  of  previous  analyses  of  DPG  methods  for  other 
equations.  Reliability  and  efficiency  estimates  for  an  error  indicator  also  follow.  Finally, 
the  equivalence  of  stability  for  various  formulations  of  the  same  Maxwell  problem  is  proved, 
including  the  strong  form,  the  ultraweak  form,  and  a  spectrum  of  forms  in  between. 


1.  Introduction 

When  a  domain  17  is  partitioned  into  elements,  a  function  in  a  Sobolev  space  like  H (curl,  17) 
or  H{ div,  17)  has  continuity  constraints  across  element  interfaces,  e.g,  the  former  has  tangen¬ 
tial  continuity,  while  the  latter  has  continuity  of  its  normal  component.  If  these  continuity 
constraints  are  removed  from  the  space,  then  we  obtain  “broken”  Sobolev  spaces.  Discontinu¬ 
ous  Petrov  Galerkin  (DPG)  methods  introduced  in  [14,  16]  used  spaces  of  such  discontinuous 
functions  in  broken  Sobolev  spaces  to  localize  certain  computations.  The  studies  in  this  paper 
begin  by  clarifying  this  process  of  breaking  Sobolev  spaces.  This  process,  sometimes  called 
hybridization,  has  been  well  studied  within  a  discrete  setting.  For  instance,  the  hybridized 
Raviart-Thomas  method  [5,  11,  28]  is  obtained  by  discretizing  a  variational  formulation  and 
then  removing  the  continuity  constraints  of  the  discrete  space,  i.e.,  by  discretizing  first  and 
then  hybridizing.  In  contrast,  in  this  paper,  we  identify  methods  obtained  by  hybridizing 
first  and  then  discretizing,  a  setting  more  natural  for  DPG  methods.  We  then  take  this  idea 
further  by  connecting  the  stability  of  formulations  with  broken  spaces  and  unbroken  spaces, 
leading  to  the  first  convergence  proof  of  a  DPG  method  for  Maxwell  equations. 

The  next  section  (Section  2)  is  devoted  to  a  study  of  the  interface  spaces  that  arise  when 
breaking  Sobolev  spaces.  These  infinite-dimensional  interface  spaces  can  be  used  to  connect 
the  broken  and  the  unbroken  spaces.  The  main  result  of  Section  2,  contained  in  Theorem  2.3, 
makes  this  connection  precise  and  provides  a  surprisingly  elementary  characterization  (by 
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duality)  of  the  natural  norms  on  these  interface  spaces.  This  theorem  can  be  viewed  as  a 
generalization  of  a  similar  result  in  [29] . 

Having  discussed  breaking  spaces,  we  proceed  to  break  variational  formulations  in  Sec¬ 
tion  3.  The  motivation  for  the  theory  in  that  section  is  that  some  variational  formulations 
set  in  broken  spaces  have  another  closely  related  variational  formulation  set  in  their  unbroken 
counterpart.  This  is  the  case  with  all  the  formulations  on  which  the  DPG  method  is  based. 
The  main  observation  of  Section  3  is  a  simple  result  (Theorem  3.1)  which  in  its  abstract 
form  seems  to  be  already  known  in  other  studies  [22],  In  the  DPG  context,  it  provides  suffi¬ 
cient  conditions  under  which  stability  of  broken  forms  follow  from  stability  of  their  unbroken 
relatives.  As  a  consequence  of  this  observation,  we  are  able  to  drastically  simplify  many  pre¬ 
vious  analyses  of  DPG  methods.  The  content  of  Sections  2  and  3  can  be  understood  without 
reference  to  the  DPG  method. 

A  quick  introduction  to  the  DPG  method  is  given  in  Section  4,  where  known  conditions 
needed  for  a  priori  and  a  posteriori  error  analysis  are  also  presented.  One  of  the  conditions  is 
the  existence  of  a  Fortin  operator.  Anticipating  the  needs  of  the  Maxwell  application,  we  then 
present,  in  Section  5,  a  sequence  of  Fortin  operators  for  H1  (K) ,  H (curl,  K)  and  H(div,  K), 
all  on  a  single  tetrahedral  mesh  element  K.  They  are  constructed  to  satisfy  certain  moment 
conditions  required  for  analysis  of  DPG  methods.  They  fit  into  a  commuting  diagram  that 
help  us  prove  the  required  norm  estimates  (see  Theorem  5.1). 

Time- harmonic  Maxwell  equations  within  a  cavity  are  considered  afterward  in  Section  6. 
Focusing  first  on  a  simple  DPG  method  for  Maxwell  equation,  called  the  primal  DPG  method, 
we  provide  a  complete  analysis  using  the  tools  developed  in  the  previous  section.  To  under¬ 
stand  one  of  the  novelties  here,  recall  that  the  wellposedness  of  the  Maxwell  equations  is 
guaranteed  as  soon  as  the  excitation  frequency  of  the  harmonic  wave  is  different  from  a 
cavity  resonance.  However,  this  wellposedness  is  not  directly  inherited  by  most  standard 
discretizations,  which  are  often  known  to  be  stable  solely  in  an  asymptotic  regime  [25].  The 
discrete  spaces  used  must  be  sufficiently  fine  before  one  can  even  guarantee  solvability  of  the 
discrete  system,  not  to  mention  error  guarantees.  Furthermore,  the  analysis  of  the  standard 
finite  element  method  does  not  clarify  how  fine  the  mesh  needs  to  be  to  ensure  that  the  stable 
regime  is  reached.  In  contrast,  the  DPG  schemes,  having  inherited  their  stability  from  the 
exact  equations,  are  stable  no  matter  how  coarse  the  mesh  is.  This  advantage  is  striking 
when  attempting  robust  adaptive  meshing  strategies. 

Mastery  of  the  proliferation  of  formulations  for  the  one  Maxwell  boundary  value  problem 
is  another  focus  of  Section  6.  One  may  decide  to  treat  individual  equations  of  the  Maxwell 
system  differently,  e.g.,  one  equation  may  be  imposed  strongly,  while  another  may  be  imposed 
weakly  via  integration  by  parts.  Mixed  methods  make  a  particular  choice,  while  primal 
methods  make  a  different  choice.  We  will  show  (see  Theorem  6.3)  that  from  the  standpoint 
of  wellposedness,  one  is  no  different  from  another,  i.e.,  if  any  one  of  the  six  formulations 
considered  in  Section  6  is  stable,  then  all  the  remaining  are  guaranteed  to  be  stable.  This 
result  is  an  interesting  application  of  the  closed  range  theorem.  However,  when  the  DPG 
methodology  is  applied  to  discretize  these  formulations,  the  numerical  results  reported  in 
Section  7,  show  that  the  various  methods  do  exhibit  differences.  This  is  because  the  functional 
settings  are  different  for  different  formulations,  i.e.,  convergence  to  the  solution  occurs  in 
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different  norms.  Section  7  also  provides  results  from  numerical  investigations  on  issues  where 
the  theory  is  currently  silent. 


2.  Breaking  Sobolev  spaces 

In  this  section,  we  discuss  precisely  what  we  mean  by  breaking  Sobolev  spaces  using  a 
mesh.  We  will  define  broken  spaces  and  interface  spaces  and  prove  a  duality  result  that 
clarifies  the  interplay  between  these  spaces.  We  work  with  infinite-dimensional,  but  rnesh- 
dependent  spaces  on  an  open  bounded  domain  1?  c  R3  with  Lipschitz  boundary.  The  mesh, 
denoted  by  Qh,  is  a  disjoint  partitioning  of  Q  into  open  elements  K  such  that  the  union  of 
their  closures  is  the  closure  of  Q.  The  collection  of  element  boundaries  dK  for  all  K  e  Qh,  is 
denoted  by  dQh-  We  assume  that  each  element  boundary  dK  is  Lipschitz.  The  shape  of  the 
elements  is  otherwise  arbitrary  for  now. 

We  focus  on  the  most  commonly  occurring  first  order  Sobolev  spaces  of  real  or  complex 
valed  functions,  namely  i71(f 7),  i7(div,  17),  and  77(curl,  17).  Their  broken  versions  are  defined, 
respectively,  by 

H\Qh)  =  {u  e  L2(! 7)  :  u\K  e  H\K),  K  e  Qh}  =  []  H\K), 

KeS 4 

H (curl,  Qh)  =  {Ee  (L2(Q))N  :  E\K  e  H{ curl,  K),  K  e  Qh)  =  ]~^  (curl,  K), 

I<eQh 

i7(div,  17/j)  =  {a  e  (L2(Q))N  :  a\K  e  H(div,K),  Ke  Qh }  =  H(div,K). 

I<eQh 

As  these  broken  spaces  contain  functions  with  no  continuity  requirements  at  element  inter¬ 
faces,  their  discretization  is  easier  than  globally  conforming  spaces. 

To  recover  the  original  Sobolev  spaces  from  these  broken  spaces,  we  need  traces  and  inter¬ 
face  variables.  First,  let  us  consider  these  traces  on  each  element  K  in  Qh- 

trgrad  u  =  u\ 8K  U  6  Hl(K), 

tr^rl  T  E  =  (■ nK  x  E)  x  nK\dK  E  e  H (curl,  K), 

tr^rlH  E  =  uk  x  E \sk  E  e  curl,  K), 

trfj)v  cr  =  a\8K  ■  nK  as  H(  div,  K). 

Here  and  throughout  uk  denotes  the  unit  outward  normal  on  dK  and  is  often  simply  written 
as  n.  Both  uk  and  these  traces  are  well  defined  almost  everywhere  on  dK ,  thanks  to  our 
assumption  that  dK  is  Lipschitz.  The  operators  trgraci,  trcuri)T,  trcuriH,  and  trdw  perform  the 
above  trace  operation  element  by  element  on  each  of  the  broken  spaces  we  defined  previously, 
thus  giving  rise  to  linear  maps 

trgrad  :  H](Qh)  ->  Y\  H1/2(dK),  trcuri>T  :  i7(curl,  Qh)  — >  Y\  H  1/2 (curl,  dK), 

KeOh  Kefih 

lLCUri,H  •  77(curl,  Qff)  *  Yl  H~1/2{div,dK),  trdiv  :  H(div,  Qh)  ]^[  H~1/2(dK). 

Keflh  KeQh 

It  is  well  known  that  these  maps  are  continuous  and  surjective.  An  element  of  WkeQ,  H~^2(dK) 
is  expressed  using  notations  like  n  ■  a  or  dn  (even  when  a  itself  has  not  been  assigned  any 
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separate  meaning)  that  are  evocative  of  their  dependence  on  the  interface  normals.  Similarly, 
the  elements  of  the  other  trace  map  codomains  are  expressed  using  notations  like 

nxE  =  E^e  ]^[  FT1/2(div,  dK),  (n  x  E)  x  n  =  ET  e  Y\  H^2  (curl,  dK). 

Keflh  KeQh 

Next,  we  need  spaces  of  interface  functions.  We  use  the  above  trace  operators  to  define 
them,  after  cautiously  noting  two  issues  that  can  arise  on  an  interface  piece  /  =  dK+  n  dK~ 
shared  by  two  mesh  elements  K-  in  17/, .  First,  functions  in  the  range  of  trgrad  when  restricted 
to  F  is  generally  multivalued,  and  we  would  like  our  interface  functions  to  be  single  valued  in 
some  sense.  Second,  the  range  of  the  remaining  trace  operators  consist  of  functionals  whose 
restrictions  to  /  are  in  general  undefined.  The  following  definitions  circumvent  these  issues. 

Hl/2(dfih)  =  trgrad  Hl(Q),  iW1/2( div,  fih)  =  trcurl!T  H( curl,  17), 

1,2 (curl,  3(2^)  =  trcuriH  H( curl,  17),  H^l^2(dfih)  =  trd;v  H(div,  17). 

Since  trgrad  Hl(fi)  C  trgrad  i71(l7/).),  one  could  endow  the  subspace  with  the  subspace  topology 
(and  similarly  for  other  spaces).  But  we  proceed  differently  and  norm  each  of  the  above 
interface  spaces  by  these  quotient  norms: 


(2.1a) 

(2.1b) 


(2.1c) 

(2-ld) 


“II  m/2{8Qh) 


I  -^t  II  H~1/2(div,8(2fl) 
E-\\\H-V2(curl,8nh) 


8n\H-y2{8Qh) 


inf 

inf 

EeH  (curl,  Q)  ntr^,^  {EH 1 


II 


\\E\ 


i/(curl,i?)  i 


inf  ,  .  ||-®||mCuri,i2)) 

EeH{Cnr\,Q)r>tr~fiXT{ET} 


inf 

ogH (div, Q)  n trjl  { an  } 


HcrlliT(div,l7)- 


These  are  indeed  quotient  norms  because  the  infimums  are  over  cosets  generated  by  kernels 
of  the  trace  maps.  E.g.,  if  w  is  any  function  in  i71(l7)  such  that  trgrad  vj  =  u,  then  the  set 
where  the  minimization  is  carried  out  in  (2.1a),  namely  i/1(l7)  n  ti'grad{h},  equals  the  coset 
w+\\KGQh  ^1( K ).  Note  that  every  element  of  this  coset  is  an  extension  of  u.  For  this  reason, 
such  norms  are  also  known  as  the  “minimum  energy  extension”  norms.  For  an  alternate  way 
to  characterize  the  interface  spaces,  see  [30]. 


Remark  2.1.  The  quotient  norm  in  (2. Id)  appeared  in  the  literature  as  early  as  [29].  The 
word  “hybrid”  that  appears  in  their  title  was  used  to  refer  to  situations  where,  to  quote  [29], 
“the  constraint  of  interelement  continuity  has  been  removed  at  the  expense  of  introducing 
a  Lagrange  multiplier.”  The  quote  also  summarizes  the  discussion  of  this  section  well.  The 
above  definitions  of  our  four  interface  spaces  are  thus  generalizations  of  a  definition  in  [29] 
and  each  can  be  interpreted  as  an  appropriate  space  of  Lagrange  multipliers. 


We  now  show  by  elementary  arguments  that  the  quotient  norms  on  the  two  pairs  of  trace 
spaces 

{H^2{dK),H-^2(dK)}  and  {FT1/2 (curl,  dK),  H^2 (div,  dK)}, 

are  dual  to  each  other.  The  duality  pairing  in  any  Hilbert  space  X,  namely  the  action  of  a 
linear  or  conjugate  linear  (antilinear)  functional  x'  e  X'  on  x  e  X  is  denoted  by  {x\x)x'xx 
and  we  omit  the  subscript  in  this  notation  when  no  confusion  can  arise.  We  also  adopt  the 
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convention  that  when  taking  supremum  over  vector  spaces  (such  as  in  the  next  result)  the 
zero  element  is  omitted  tacitly. 

Lemma  2.2.  The  following  identities  hold  for  any  n-a  in  H~1/l2(dK)  and  any  u  in  H1/l2( dK ). 


(2.2a) 

(2.2b) 


II  II  #-1/2(010  ~  SUP 

ugH 1(K) 


I  (Vn,u) 


=  sup 


\(Vn,U)  | 


\\u\\  Hl{K)  ueW/^dK)  INI  m/2(dK) 

\(n-(T,u)\  |  <dn ,  u)  | 


\\u\\m/H8K)=  SUP  ' '  I,  I, '  ' 1  =  sup 

<reH(div,K)  \\a\\H(dn,K)  &neH-V2(dK)  Wari\\  H~V2(dK) 


The  next  two  identities  hold  for  any  EH  in  H  1/2(div,d/i)  and  any  FT  in  H  ^(curl,  dK). 


(2.2c) 


WE^Wn-^tdiv^K)  ~  SUP 

FeH(curl,K) 


=  sup 
i'Ye.ff-1/2  (curl, SIC) 


'  \\H(cut\,K) 

\(E^Ft) 


II F 


T II  H~1t2(curl,8K) 


(2. 2d) 


\\Et\\h-i/2(cu1.\  8k)  -  sup 

EeH(curl,K) 


sup 


|<n  x  E,  FT 


\\E\\h  (curl, K) 

\<E^Fr)\ 


E-ieH-iPidivJK) 


II E 


'-\V\H~W2  (div,dK) 


Proof.  The  Hrst  identity  is  proved  using  an  equivalence  between  a  Dirichlet  and  a  Neu¬ 
mann  problem.  The  Dirichlet  problem  is  the  problem  of  finding  a  £  H(div,  K),  given 
an  £  H^1^2(dK),  such  that 


(2.3) 


{n  ■  a  =  n  ■  a, 

—  grad(div  a)  +  cr  =  0, 


The  Neumann  problem  finds  w  £  Hl(K)  satisfying 

{d iv 

=  CJn’ 

— div(gradtc)  +  w  =  0, 


on  dK, 
in  K. 


on  dK, 
in  K. 


It  is  immediate  that  problems  (2.4)  and  (2.3)  are  equivalent  in  the  sense  that  w  solves  (2.4) 
if  and  only  if  a  =  gradtc  solves  (2.3)  and  moreover  ||wi||/o(/o  =  ||o'||ij(div,A')-  It  is  also  obvious 
from  the  calculus  of  variations  that  among  all  H{ div,  It ) -extensions  of  dn,  the  solution  of  (2.3) 
has  the  minimal  H(div,/i)  norm  (i.e. ,  a  is  the  “minimum  energy  extension”  referred  to 
earlier),  so 


PnllH^tdK)  ~  llCJII.H'(div,10  ~  WIh1^) 

\(dn,v)\ 

=  sup  — - , 

vem(K)  IMh1  (K) 


sup 

vgH 1(K) 


|  (grad  w,  grad  v) k  +  ( w,v)k 
MIhPk) 
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where  we  used  the  variational  form  of  (2.4)  in  the  last  step.  (Here  and  throughout,  we  use 
(•,  -)k  to  denote  the  inner  product  in  L2(K)  or  its  Cartesian  products.)  This  proves  the  first 
equality  of  (2.2a). 

Next,  analogous  to  (2.3)  and  (2.4),  we  set  up  another  pair  of  Dirichlet  and  Neumann 
problems.  The  first  problem  is  to  find  u  in  Hl(K ),  given  any  u  £  ELl/2(dK),  such  that 


(2.5) 


u  =  u,  on  dK , 

— div(gradu)  +  u  =  0,  in  K. 


The  second  is  to  find  r  in  H( div,  K )  such  that 


(2.6) 


div  t  =  u,  on  dK, 

—  grad(div  r)  +  r  =  0,  in  K. 


The  solution  u  of  (2.5)  has  the  minimal  H 1(K)  norm  among  all  extensions  of  u  into  Hl{K), 
i.e.,  \\u\\hi/2(8K)  =  \\u\\h1(k)-  Thus  (an,  u)/\\u\\m/*(dK)  =  (?n,  u)/\\u\\hi(k),  so  taking  the 
supremum  over  all  u  in  i71//2(d.A),  we  obtain 


sup 

ueH1/2(8K) 


\(Vn,u)\ 

l|u||_ffi/2(dx) 


^  sup 

ueH1(K) 


\(Vn,u)\ 

I \u\\m(K) 


Since  the  reverse  inequality  is  obvious  from  the  definition  of  the  quotient  norm  in  the  denom¬ 
inator,  we  have  established  the  second  identity  of  (2.2a).  To  prove  (2.2a),  we  begin,  as  above, 
by  observing  that  r  is  the  solution  to  the  Neumann  problem  (2.6)  if  and  only  if  u  =  divr 
solves  the  Dirichlet  problem  (2.5).  Moreover,  ||r||//(diV]x)  =  II u Hi?1  (K)-  Hence 


^\\mi2{dK) 


u\\  H'iK) 


|(div  r,  div  p)k  +  (t,  p)k 

r\\H(div,K)  =  ST-lP  - j]— j] - 

peH  (div, K)  llP!lif(div,A:) 


sup 

peH(div,K) 


\(n  ■  P,  u) | 

\p\\H(div,K) 


where  we  have  used  the  variational  form  of  (2.6)  in  the  last  step.  The  proof  of  (2.2b)  can 
now  be  completed  as  before. 

We  follow  exactly  the  same  reasoning  for  the  H{ curl)  case,  summarized  as  follows:  On  one 
hand,  the  norm  of  an  interface  function  equals  the  norm  of  a  minimum  energy  extension,  while 
on  the  other  hand,  it  equals  the  norm  of  the  inverse  of  a  Riesz  map  applied  to  a  functional 
generated  by  the  interface  function.  The  minimum  energy  extension  that  yields  the  interface 
norm  ||AH||^_i/2(div  is  now  the  solution  of  the  Dirichlet  problem  of  finding  E  e  H{ curl,  K ) 
satisfying 


(2.7) 


n  x  E  =  E_ |,  on  dK, 

curl(curl  E)  +  E  =  0,  in  K, 


while  the  inverse  of  the  Riesz  map  applied  to  the  functional  generated  by  is  obtained  by 
solving  the  Neumann  problem 


(2.8) 


n  x  (curlF)  =  AH, 
curl  (curl  F )  +  F  =  0, 


on  dK 
in  K. 
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Again,  the  two  problems  are  equivalent  in  the  sense  that  F  solves  (2.8)  if  and  only  if  E  = 
curlF  solves  (2.7).  Moreover,  \\E\\H(cuiiK)  =  \\F\\ h(cuiI,k)-  Hence 

I -S'h  II (div,5/4T)  =  W^\\h(cuiI,K)  =  II F\\ //(curl, K) 

|  (curl  F,  curl  G)k  +  (F,  G)k\ 

sup  - \FF\ - 

GeH(cm\,K)  IMIir(curl,.fO 


sup 

GeH(cur\,K) 


\<E^G>  | 

II G II  //(curl,  K) 


The  proof  of  (2.2c)  follows  from  this.  The  proof  of  (2. 2d)  is  similar  and  is  left  to  the 
reader.  □ 


Let  us  return  to  the  product  spaces  like  L71(l7/l),  H( curl,  f2ff),  and  Lf(div,  f?h)-  Any  Hilbert 
space  V  that  is  the  Cartesian  product  of  various  Hilbert  spaces  V(K)  is  normed  in  the 
standard  fashion, 

V=  El  VW’  Hlv  =  S  H^ll v(Ky 

Kefih  Kefih 

where  vk  denotes  the  il -component  of  any  v  in  V .  The  dual  space  V'  is  the  Cartesian  product 
of  component  duals  V(K)'.  Writing  an  £  e  V'  as  t(v)  =  YjKsQh  &k(vk),  where  Ik  £  V(K)', 
it  is  elementary  to  prove  that  ||£|| y,  =  2 KeQh  ll^lly(K)'>  ue., 

(2.9)  LH  y  (  sup  !M  V 

\veV  MW  J  KeQh  \vkeV(K)  \\Vk\\v(K)  J 

Some  of  our  interface  spaces  have  such  functionals,  e.g.,  the  function  an  in  iL~1//2(dL?/j)  gives 
rise  to  t(v)  =  (a n,v}h  where 

(pn,V)h=  2  (VmV)H-i/2(8K)xHl/2(8K), 

KeS7h 

is  a  functional  acting  on  v  s  iL1(l Iff)  which  is  the  sum  of  component  functionals  £k(v)  = 
(<rn,  v} H -i/2 r 8 k)xH1/2( 3 K)  acting  on  vk  =  v\k  over  every  K  e  Other  functionals  like 
(fE_ | ,  F)8  are  defined  similarly.  We  are  now  ready  to  state  a  few  basic  relationships  between 
the  interface  and  broken  spaces. 


Theorem  2.3.  The  following  identities  hold  for  any  interface  space  function  dn  in  H  1/2(dl7/l), 
u  in  H1^2(df2h),  E^  in  it~1//2(div,  dK),  and  FT  in  iL-1/2(curl,  dK). 


(2.10a) 

(2.10b) 

(2.10c) 

(2.10d) 


llcrn||_ff-l/2(^)  —  SUP 


\(Vn,u)h\ 


\\u\\m/^8Qh)  =  SUP 


uemlQb) 

\{n  •  a,  u)h\ 


WE^WH-^tdiVtdQh)  -  sup 

FeH(curl,ah) 


a<=H(div,f2h)  ll^llif  (div,,!?^) 

\<E^,F)h  | 


II A I 


if  (curl, 1?^) 


l-*'V||.ff-i/2(curi  d1f4)  -  sup 

EeH{c  url,fih) 


\{n  x  E,FT)h | 


\\E\ 


if  (curl,  4?^,) 
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For  any  broken  space  function  v  e  H1(f2f1),r  £  //( div,  17/,),  and  F  e  //(curl,  17/,), 


(2.11a) 

veH1{Q) 

(Ani 

O 

II 

S> 

V&ne  H~V2(dnh), 

(2.11b) 

t  £  //(div,  17) 

(t  • 

n,  u)h  =  0 

Vue  //1/2(dl4), 

(2.11c) 

F  £  //(curl,  n) 

*=*  <fh 

)  F)h  =  0 

V/h£  //~1/2(div,  17/,) 

Proof.  The  identities  immediately  follow  from  Lemma  2.2  and  (2.9).  The  proofs  of  the  three 
equivalences  in  (2.11)  are  similar,  so  we  will  only  detail  the  last  one.  If  F  is  in  //(curl,  17), 
then  choosing  any  E  £  //(curl,  17)  such  that  trcur]T  E  =  FH  and  integrating  by  parts  over 
entire  17, 


(curl  E,F)q  +  (E,cut1F)q  =  0 

because  of  the  boundary  conditions  on  F  on  dft.  Now,  if  the  left  hand  side  is  integrated  by 
parts  again,  this  time  element  by  element,  then  we  find  that  (E_^,  F)/,  =  0. 

Conversely,  given  that  (F_,,  F)/,  =  0  for  any  F  in  //( curl,  17/,),  consider  curl  F  e  (X^l?)3)'. 
As  a  distribution,  curlF  acts  on  <f>  e  D(l7)3,  and  satisfies 

(curl  F)(4>)  =  (F,  curl  </>)/?  =  (curlF,  </>)/,  -  (n  x  <f>,  F)h  =  (curl  F,<j>)h, 

where  we  have  integrated  by  parts  element  by  element  and  denoted 

(’’  ')h  =  Yj  (’> ') K ■ 

KeS lh 

This  notation  also  serves  to  emphasize  that  the  term  curl  F  appearing  on  the  right-hand  side 
above  is  a  derivative  taken  piecewise,  element  by  element.  Clearly  curlF^-  is  in  L2(K )3  for  all 
K  £  17/,  since  F  £  //(curl,  17/,) ,  so  the  distribution  curl/1  is  in  L2(l 7).  Having  established  that 
F  £  //(curl,  17),  we  may  now  integrate  by  parts  to  get  (n  x  E,  F>#m/2(diV)5J7\x#m/2(curl  = 

(curl  E,F)q  +  (F,  curlF)^  =  (n  x  F,  F)/,  =  0  for  all  E  £  //(curl,  17).  This  shows  that  the 
trace  (n  x  F)  x  n\gQ  =  0,  i.e,  F  £  //(curl,  17).  □ 

Remark  2.4.  While  nMl  and  nAs«,  H'/\dK)  are  dual  to  each  other,  our 

interface  spaces  /7~1//2(dl7/,)  and  /71//2(dl7/,)  are  not  dual  to  each  other  in  general. 

Remark  2.5.  Equivalences  analogous  to  (2.11)  hold  with  interface  subspaces 

H1/2{dQh)  =  trgrad  H~1,2( div,  Qh)  =  trcurlH  //(curl,  C), 

//~l/2(curl,  dfih)  =  trcuri;T  //(curl,  17),  H~]'2(dflh)  =  trdiv  //(div,  17). 

By  a  minor  modification  of  the  arguments  in  the  proof  in  Theorem  2.3,  we  can  prove  that 
for  any  v  £  //1(l7/l),  r  £  //(div,  17/,.),  and  F  £  //(curl,  17/,), 


(2.12a) 

u  £  Hl{Q)  <= = =>  (an,v)h  =  0 

V^E//-1/2^/,), 

(2.12b) 

t  e  H (div,  17)  <^=>  (r  •  n,  u)/,  =  0 

Vii  £  //1/2(dl7/, ) , 

(2.12c) 

F  £  //(curl,  17)  (E^,F)h  =  0 

VFH  £  Z/_1//2(div,  17/, 
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3.  Breaking  variational  forms 


The  goal  in  this  section  is  to  investigate  in  what  sense  a  variational  formulation  can  be 
reformulated  using  broken  spaces  without  losing  stability.  We  will  describe  the  main  result 
in  an  abstract  setting  first  and  close  the  section  with  simple  examples  that  use  the  results  of 
the  previous  section. 

Let  Xo  and  Y  denote  two  Hilbert  spaces  and  let  Yq  be  a  closed  subspace  of  Y.  For 
definiteness,  we  assume  that  all  our  spaces  in  this  section  are  over  C  (but  our  results  hold  also 
for  spaces  over  M).  In  the  examples  we  have  in  mind,  Y  will  be  a  broken  space,  while  Yq  will  be 
its  unbroken  analogue  (but  no  such  assumption  is  needed  to  understand  the  upcoming  results 
abstractly).  The  abstract  setting  involves  a  continuous  sesquilinear  form  bo  :  Xq  x  Y  — *  C 
satisfying  the  following  assumption. 


Assumption  1.  There  is  a  positive  constant  cq  such  that 

\bo(x,y)\ 


cofHxo  sS  sup 
yeY0 


\\y\\Y 


MxeX0. 


It  is  a  well-known  result  of  the  Babuska  and  Necas  [1,  26]  that  Assumption  1  together  with 
triviality  of 


(3.1)  Zo  =  {y  e  Yo  :  bo(x,  y)  =  0  for  all  x  £  Ao} 

guarantees  wellposedness  of  the  following  variational  problem:  Given  £  e  Yq  (the  space  of 
conjugate  linear  functionals  on  Yo),  find  x  £  Xo  satisfying 

(3.2)  bo(x,  y)  =  £{y)  Vy  e  Y0. 


When  Zq  is  non-trivial,  we  can  still  obtain  existence  of  a  solution  x  provided  the  load  func¬ 
tional  £  satisfies  the  compatibility  condition  £(z)  =  0  for  all  z  £  Zq-  In  (3.2),  the  trial  space 
Xo  need  not  be  the  same  as  the  test  space  Yq. 

To  describe  a  “broken”  version  of  (3.2),  we  need  another  Hilbert  space  X,  together  with  a 
continuous  sesquilinear  form  b  :  X  x  Y  — >  C.  In  applications  Y  and  X  will  usually  be  set  to 
a  broken  Sobolev  space  and  an  interface  space,  respectively.  Define 

b(  (x,  x),y)  =  b0(x,  y)  +  b(x,  y). 


Clearly  b  :  X  x  Y  — *  C  is  continuous  where 

(3.3)  X  =  X0  x  X 

is  a  Hilbert  space  under  the  Cartesian  product  norm.  Now  consider  the  following  new  broken 
variational  formulation:  Given  l  e  Yf,  find  x  £  Xq  and  x  £  X  satisfying 

(3.4)  b{(x,x),  y)  =  £(y)  Vy  e  Y. 

The  close  relationship  between  problems  (3.4)  and  (3.2)  is  readily  revealed  under  the  following 
assumption. 


Assumption  2.  The  spaces  Yo,Y,  and  X  satisfy 

(3.5)  Yq  =  {y  e  Y  :  b(x,  y)  =  0  for  all  x  £  X] 
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and  there  is  a  positive  constant  c  such  that 

-||*n  ^  IKK  y)\  - 

(3.6)  c  s  U  <  sup-r-i -  Visa. 

X  yeY  My 

Under  this  assumption,  we  present  a  simple  result  which  shows  that  the  broken  form  (3.4) 
inherits  stability  from  the  original  unbroken  form  (3.2).  A  very  similar  such  abstract  result 
was  formulated  and  proved  in  [22,  Appendix  A]  and  used  for  other  applications.  Our  proof 
is  simple,  unsurprising,  and  uses  the  same  type  of  arguments  from  the  early  days  of  mixed 
methods  [5,  p.  40]:  stability  of  a  larger  system  can  be  obtained  in  a  triangular  fashion  by  first 
restricting  to  a  smaller  subspace  and  obtaining  stability  there,  followed  by  a  backsubstitution- 
like  step. 


Theorem  3.1.  Assumptions  1  and  2  imply 


ci  ||  (z,®)  || x  sS  sup 


IK  (x,x),  y)  | 


yeY  \\y\\Y 


where  c\  is  defined  by 


1  1,1  /||60 


„2  „2  1  pa 

1  c0  ' 


Co 


+ 1 


Moreover,  if  Z  =  {y  e  Y  :  b(  (. x ,  x),  y)  =  0  for  all  x  e  Xq  and  x  in  X},  then 


Z  =  Z0. 


Consequently,  if  Zq  =  {0},  then  (3.4)  is  uniquely  solvable  and  moreover  the  solution  compo¬ 
nent  x  from  (3.4)  coincides  with  the  solution  of  (3.2). 


Proof.  We  need  to  bound  ||x||x0  and  ||Klx-  First, 


cofUxq  sup 
yo  sY0 

^  sup 
yoeYo 


\bp(x,y)\ 

\\y\\  y 

I  bp{x,y)  +  b(x,y)\ 

\\y\\Y 


„„„  \K(x,x),  y) I 

□  Up  ..  .. 

yeY  ||y||y 


Next,  to  bound  ||x||y>  using  (3.6)  of  Assumption  2, 


by  Assumption  1, 


by  Assumption  2,  (3.5) 


as  Yq  c  Y. 


C  ||x|| g  ^  sup 
yeY 


\b{x,y)\ 

hh 


_  IK  (x,x),  y)  -  b0(x,y)\ 
sup  ..  || 

yeY  \\y\\Y 


<  INI  \\x\\x0  +  sup 

yeY 


\b((x,x),  y)  | 


\\y\\Y 


where  |NI  is  the  smallest  number  C  for  which  the  inequality  \bo(x,y)\  ^  C,||®||x0||y||Y  holds 
for  all  x  e  Xq  and  all  yeY.  Using  the  already  proved  bound  for  ||x||x0  hi  the  last  inequality 
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and  combining, 


\(x,x)\\x  —  \\x\\x0  + 


"X 


\Cq  yeY  My 


from  which  the  inequality  of  the  theorem  follows. 
Finally,  to  prove  that  Z  =  Zq,  using  (3.5), 


y  e  Z  <^=>  bo(x,  y)  =  0  for  all  x  e  Xq  and  b(x,  y)  =  0  for  all  x  e  X 
<==>  bo(x,  y)  =  0  for  all  x  e  Xq  and  y  e  Yq, 


which  holds  if  and  only  if  y  £  Zq. 


□ 


Remark  3.2.  Note  that  in  the  proof  of  the  inf-sup  condition,  we  did  not  fully  use  (3.5).  We 
only  needed  Yq  c  {y  e  Y  :  b(x,y )  =  0  for  all  x  e  X}.  The  reverse  inclusion  was  needed  to 
conclude  that  Z  =  Zq. 


Remark  3.3.  It  is  natural  to  ask,  in  the  same  spirit  as  Theorem  3.1,  if  the  numerical  solutions 
of  DPG  methods  using  discretizations  of  the  broken  formulations  coincide  with  those  of 
discretizations  of  the  original  unbroken  formulation.  A  result  addressing  this  question  is 
given  in  [3,  Theorem  2.6]. 

In  the  remainder  of  this  section,  we  illustrate  how  to  apply  this  theorem  on  some  examples. 

Example  3.4  (Primal  DPG  formulation).  Suppose  /  e  L2(I2)  and  u  satisfies 

(3.7a)  —A  u  =  f  in  17, 

(3.7b)  u  =  0  on  dfi. 

The  standard  variational  formulation  for  this  problem,  finds  u  in  i71(I?)  such  that 

(3.8)  (grad  u,  grad  v)q  =  (f,v)n  Vv  £  if1  (12). 

This  form  is  obtained  by  multiplying  (3.7a)  by  v  £  i71(I2)  and  integrating  by  parts  over  the 
entire  domain  12.  If  on  the  other  hand,  we  multiply  (3.7a)  byare  F71(I2/l)  and  integrate  by 
parts  element  by  element,  then  we  obtain  another  variational  formulation  proposed  in  [17]: 
Solve  for  u  in  i71(I2)  as  well  as  a  separate  unknown  an  e  H~1^2(df2h)  (representing  the 
fluxes  n  ■  grad  u  along  mesh  interfaces)  satisfying 

(3.9)  (grad it,  grad y)h  +  (&n,y)h  =  {f,y)n  Vy  e  Hl{flh). 

We  can  view  this  as  the  broken  version  of  (3.8)  by  setting 

X0  =  H1(fl),  Yq  =  H\Q ), 

X  =  H-y\dS2h),  Y  =  H1(nh), 

bo(u,y)  =  (grad  rq  grad  y)h,  b(an,y)  =  < an,y)h . 
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Figure  1 .  Chains  of  implications  of  inf-sup  conditions 

For  these  settings,  the  conditions  required  to  apply  Theorem  3.1  are  verified  as  follows. 

Coercivity  of  bo(-,  •)  on  To  =>  Assumption  1  holds. 

Theorem  2.3,  (2.10a)  =^>  (3.6)  of  Assumption  2  holds  with  c  =  1. 

Theorem  2.3,  (2.11a)  =>  (3.5)  of  Assumption  2  holds. 

Noting  that  Zq  =  {0},  an  application  of  Theorem  3.1  implies  that  problem  (3.9)  is  wellposed. 
This  wellposedness  result  also  shows  that  (3.9)  is  uniquely  solvable  with  a  more  general 
right-hand  side  /  in  H1  (12^)'. 

An  alternate  (and  longer)  proof  of  this  wellposedness  result  can  be  found  in  [17].  The 
classical  work  of  [29]  also  uses  the  spaces  Hl(flh)  and  17-1/2(dl?fc),  but  proceeds  to  develop  a 
Bubnov-Galerkin  hybrid  formulation  different  from  the  Petrov-Galerkin  formulation  (3.9).  /// 

Example  3.5  (Many  formulations  of  an  elliptic  problem).  Considering  a  model  problem  in¬ 
volving  diffusion,  convection,  and  reaction  terms,  we  now  show  how  to  analyze,  all  at  once, 
its  various  variational  formulations.  The  diffusion  coefficient  a  =  a-1  :  12  — *  M3x3  is  a  sym¬ 
metric  matrix  function  which  is  uniformly  bounded  and  positive  definite  on  12,  the  convection 
coefficient  is  /3  e  L°°(12)3  which  satisfies  div(a/3)  =  0,  and  reaction  is  incorporated  through  a 
non- negative  7  e  L°°(12).  The  classical  form  of  the  equations  on  12  are  a  =  a  grad  u  +  aj3u  +  f± 
and  —  div  a  +  'yu  =  (for  some  given  f\  e  L2(12)3  and  f2  e  L2(12))  together  with  the  bound¬ 
ary  condition  u\sq  =  0.  This  can  be  written  in  operator  form  using 


A 

a 

aa  —  grad  u  —  f3u 

,  A* 

a 

aa  —  grad  u 

u 

div  a  —  'yu 

u 

div  a  —  (3  ■  a  —  'yu 

We  begin  with  the  formulation  closest  to  the  classical  form. 

Strong  form:  Let  x  =  (a,  u )  be  a  group  variable.  Set  spaces  by 

Xq  =  H{ div,  12)  x  iL1(12),  Y  =  Y0  =  L2(12)3  x  L2(12), 

and  consider  the  problem  of  finding  x  e  Aq,  given  /  e  T,  satisfying  Ax  =  f.  We  can 
trivially  fit  this  into  our  variational  framework  (3.2)  by  setting  60  to 

bo(x,y)  =  {Ax,y)n. 

Unlike  the  remaining  formulations  below,  there  is  no  need  to  discuss  a  broken  version  of 
the  above  strong  form  as  the  test  space  already  admits  discontinuous  functions.  The  next 
formulation  is  often  derived  directly  from  a  second  order  equation  obtained  by  eliminating  a 
from  the  strong  form. 
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Primal  form:  First,  set  spaces  by 

x0  =  h\s2),  y0  =  h\q), 

X  =  H-y2{dnh),  Y  =  H\Qh). 

Then,  with  b  set  to  bp(an ,  v)  =  (n  ■  a,  y)h  and  bo  set  to 

bp(u,v)  =  (a  grad u,  grad  v) h  ~  (a(3u,  gr&dv)h  +  (ju,v)a, 

the  standard  primal  formulation  is  (3.2)  and  its  broken  version  is  (3.4). 

Next,  consider  the  formulation  derived  by  multiplying  each  equation  in  the  strong  form  by  a 
test  function  and  integrating  both  equations  by  parts,  i.e.,  both  equations  are  imposed  weakly. 
It  was  previously  studied  in  [9,  15],  but  we  can  now  simplify  its  analysis  considerably  using 
Theorem  3.1. 

Ultraweak  form:  Set  group  variables  x  =  ( a,u ),  y  =  (r,  v),  x  =  ( an,u ),  and 

X0  =  L2(f2)3  x  L2(I2),  Y0  =  H{ div,  Q)  x 

X  =  H1/2(dnh)  x  H-WtfQh),  Y  =  if  (div,  Qh)  x  Hl{Qh), 
bo{x ,  y)  =  (x,  A*y)Qh ,  bL\x ,  y)  =  <<r„,  v)h  +(u,n  •  r)h. 

Formulations  (3.2)  and  (3.4)  with  bo  =  bp  and  b  =  bu  are  of  the  ultraweak  type. 

The  fourth  formulation,  well-known  as  the  mixed  form  [5],  is  derived  by  weakly  imposing  (via 
integration  by  parts)  the  first  equation  of  the  strong  form,  but  strongly  imposing  the  second 
equation. 

Dual  Mixed  form:  Set  the  spaces  by 

X0  =  if  (div,  Q)  x  L2(I2),  Yq  =  H  (div,  12)  x  L2(I2), 

X  =  H1/2(df2h),  Y  =  if  (div,  Qh)  x  L2{Q). 

The  well-known  mixed  formulation  is  then  (3.2)  with  bo  =  bp , 

bp((a,u),  (t,v)  )  =  ( aa,r)h  +  («,divr)ft  -  (/ 3u,r)h 
+  (div  a ,  v)h  -  (7 u,v)h. 

Its  broken  version  is  (3.4)  with  b  set  to  bD(u,  r)  =  (u,  n  ■  r)^. 

Note  that  the  well-known  discrete  hybrid  mixed  method  [5,  11]  is  also  derived  from  bp .  That 
method  however  works  with  a  Bubnov-Galerkin  formulation  obtained  by  breaking  both  the 
trial  and  the  test  if  (div,  12)  components,  while  above  we  have  broken  only  the  test  space.  The 
last  formulation  in  this  example  reverses  the  roles  by  weakly  imposing  the  second  equation 
of  the  strong  form  and  strongly  imposing  the  first  equation: 

Mixed  form:  Set 

X0  =  L2(I2)3  x  H1^), 

X  =  H-^dQh), 


Y0  =  L2(i?)3  x  H'({2), 
Y  =  L2(I2)3  x 
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The  dual  mixed  formulation  is  (3.2)  with  bo  =  b(f , 

b™ ( (a,  u) ,  (r,  v) )  =  (aa,  r)  q  ~  (grad  u,t)q  -  {fdu,r)n 
-  (a,  gradu)j7  -  (7 u,v)q. 

and  its  broken  version  is  (3.4)  with  b  set  to  bM(an,v)  =  (an,v)y. 

The  variational  problem  (3.2)  with  b{f  is  sometimes  called  [4]  the  primal  mixed  form  to 
differentiate  it  with  the  dual  mixed  form  given  by  5^.  The  broken  formulation  (3.4)  with 
by1  and  bM  was  called  the  mild  weak  DPG  formulation  in  [6].  Their  analysis  can  also  be 
simplified  now  using  Theorem  3.1. 

In  order  to  apply  Theorem  3.1  to  all  these  formulations,  we  need  to  verify  Assumption  1. 
This  can  be  done  for  all  the  formulations  at  once,  because  the  six  implications  displayed  in 
Figure  1  are  proved  in  [12]  for  the  model  problem  of  this  example.  We  will  not  detail  this 
proof  here  because  we  provide  full  proofs  of  similar  implications  for  Maxwell  equations  in 
Section  6  (and  this  example  is  simpler  than  the  Maxwell  case).  To  apply  these  implications 
for  the  current  example,  we  pick  a  formulation  for  which  Assumption  1  is  easy  to  prove: 
That  the  primal  form  is  coercive  b^(u,u)  ^  C\\ follows  immediately  by  integration 
by  parts  and  the  Poincare  inequality  (under  the  simplifying  assumptions  we  placed  on  the 
coefficients).  This  verifies  Assumption  1  for  the  primal  form,  which  in  turn  verifies  it  for 
all  the  formulations  by  the  above  chain  of  equivalences.  Assumption  2  can  be  immediately 
verified  for  all  the  formulations  using  either  (2.11)  or  (2.12).  Together  with  the  easily  verified 
triviality  of  Zq  in  each  case,  we  have  proven  the  wellposedness  of  all  the  formulations  above, 
including  the  broken  ones.  /// 


4.  The  DPG  method 

In  this  section,  we  quickly  introduce  the  DPG  method,  indicate  why  the  broken  spaces  are 
needed  for  practical  reasons  within  the  DPG  method,  and  recall  known  abstract  conditions 
under  which  an  error  analysis  can  be  conducted. 

Let  X  and  Y  be  Hilbert  spaces  and  let  b  :  X  x  Y  — *  C  be  a  continuous  sesquilinear  form. 
In  the  applications  we  have  in  mind,  X  will  always  be  of  the  form  (3.3)  (but  we  need  not 
assume  it  for  the  theory  in  this  section).  The  variational  problem  is  to  find  x  in  X,  given 
£  e  Y',  satisfying 

(4.1)  b{x,y)=£{y)  MyeY. 

The  DPG  method  uses  finite-dimensional  subspaces  Xy  c  X  and  Yy  c;  Y .  The  test  space 
used  in  the  method  is  a  subspace  YJ°pt  ^  Yy  of  approximately  optimal  test  functions  computed 
for  any  arbitrarily  given  trial  space  Xy.  It  is  defined  by  Y^pt  =  Ty(Xy)  where  Ty  :  Xy  —*  Yy 
is  given  by 

(4.2)  (Thz,y)Y  =  b(z,y)  VyeYy. 

Here  (•,  -)y  is  the  inner  product  in  Y,  hence  by  Riesz  representation  theorem  on  Yy,  the 
operator  Ty  is  well  defined.  The  discrete  problem  posed  by  the  DPG  discretization  is  to  find 
xy  e  Xy  satisfying 

(4.3) 


b(xh,y)=£(y)  Vy6T°pt. 
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For  practical  implementation  purposes,  it  is  important  to  note  that  7),  can  be  easily  and 
inexpensively  computed  via  (4.2)  provided  the  space  Yh  is  a  subspace  of  a  broken  space. 
Then  (4.2)  becomes  a  series  of  small  decoupled  problems  on  each  element. 

For  a  posteriori  error  estimation,  we  use  an  estimator  fj  that  actually  works  for  any  x h  in 
Xh,  computed  as  follows.  (Note  that  Xh  need  not  equal  the  solution  Xh  of  (4.3).)  First  we 
solve  for  Eh  in  1),  by 

(4.4)  (£h,y)Y  =  i(y)  -  b(xh,y),  TyeYh. 

Again,  this  amounts  to  a  local  computation  if  Yh  is  a  subspace  of  a  broken  space.  Then,  set 

fj  =  ¥h\\v- 

When  Yh.  is  a  broken  space,  the  element-wise  norms  of  Eh  serve  as  good  error  estimators  [18]. 
The  notations  y  and  Eh  (without  tilde)  refer  to  similarly  computed  quantities  with  Xh  in  place 
of  Xh-  An  analysis  of  errors  and  error  estimators  of  the  DPG  method  can  be  conducted  using 
the  following  assumption  introduced  in  [23].  In  accordance  with  the  traditions  in  the  theory 
mixed  methods  [5],  we  will  call  the  operator  IJ  in  the  assumption  a  Fortin  operator. 

Assumption  3.  There  is  a  continuous  linear  operator  17  :  Y  — ►  Yh  such  that  for  all  Wh  £  Xh 
and  all  v  e  Y, 

b(wh,  v  —  IIv )  =  0. 

Theorem  4.1.  Suppose  Assumption  3  holds.  Assume  also  that  there  is  a  positive  constant 
c\  such  that 

(4.5)  ci||a; ||  x  <  sup  - ^ -  V.t  e  X , 

yeY  MY 

and  the  set  Z  =  {y  eY  :  b(x,  y)  =  (1  for  all  x  e  X}  equals  {0}.  Then  the  DPG  method  (4.3) 
is  uniquely  solvable  for  Xh  and  the  a  priori  error  estimate 

(4.6)  || a;  —  Xh\\x  ^  — -  inf  ||x  —  Zh\\x  (quasi- optimality) 

Cl  ZfiEXfo 

holds,  where  x  is  the  unique  exact  solution  of  (4.1).  Moreover,  we  have  the  following  inequal¬ 
ities  for  any  Xh  in  Xh  and  its  corresponding  error  estimator  fj,  with  the  data- approximation 
error  osc(£)  =  ||£  o  (1  —  17)||y/. 

(4.7a)  ci|ja;  —  Xh\\x  ^  \\n\\  fj  +  osc(£), 

(4.7b)  fj  ^\\b\\\\x- xh\\x, 

(4.7c)  osc(£)  ^  ||6||  ||  1  —  I7||  min  ||.t  —  Zh\\x 

Zh.£Xh 

Here  ||il||  and  ||6||  are  any  constants  that  satisfy  ||i7y||y  ^  ||i7||||y||y  and  \b(w,y)\  ^ 
ll&lllklUI|y||y,  respectively,  for  all  w  e  X  and  yeY.  To  apply  the  theorem  to  specific 
examples  of  DPG  methods,  we  must  verify  (4.5).  This  will  usually  be  done  by  appealing  to 
Theorem  3.1  and  verifying  Assumptions  1  and  2.  The  previous  sections  provided  tools  for 
verifying  Assumptions  1  and  2.  In  the  next  section,  we  will  provide  some  tools  to  verify  the 
remaining  major  condition  in  the  theorem,  namely  Assumption  3. 


(reliability) 
( efficiency) 
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Remark  4.2.  A  proof  of  Theorem  4.1  is  available  in  existing  literature.  The  a  priori  error 
bound  (4.6)  was  proved  in  [23].  The  inequalities  of  (4.7),  useful  for  a  posteriori  error  esti¬ 
mation,  were  proved  in  [7].  In  particular,  a  reliability  estimate  slightly  different  from  (4.7a) 
(with  worse  constants)  was  proved  in  [7],  but  the  same  ideas  yield  (4.7a)  easily  (for  example, 
cf.  [8,  proof  of  Lemma  3.6]). 


Remark  4.3.  The  operator  7),  is  an  approximation  to  an  idealized  trial-to-test  operator  T  : 
X  — >  Y  given  by 

(Tx,y)Y  =  b(x,y)  MyeY. 

If  B  :  X  — >  Y'  is  the  operator  defined  by  the  form  satisfying  ( Bx)(y )  =  b(x,y )  for  all  x  e  X 
and  y  e  Y,  then  clearly  T  =  RylB,  where  Ry  :  Y  — *•  Y'  is  the  Riesz  map  defined  by 
( Ryy)(v )  =  (■ y ,  v)y  ■  In  some  examples  [15],  it  is  possible  to  analytically  compute  T  and  then 
one  may  substitute  Yj°pt  with  the  exactly  optimal  test  space  Yopt  =  T(Xh). 

Remark  4.4.  The  above-mentioned  trial-to-test  operator  T  =  R^B  should  not  be  confused 
with  another  trial-to-test  operator  S  =  ( B')~1Rx  of  [2]  (also  cf.  [20]): 

X  — Y’ 


Rx 


Ry 


X' 


B' 


Y 


Application  of  S  requires  the  inversion  of  the  dual  operator  B' . 


5.  Fortin  operators 

The  Fortin  operator  77  appearing  in  Assumption  3  is  problem  specific  since  it  depends 
on  the  form  b  and  the  spaces.  However,  there  are  a  few  Fortin  operators  that  have  proved 
widely  useful  for  analyzing  DPG  methods,  including  one  for  Y  =  771(l7 h)  and  another  for 
Y  =  77(div,  f2h),  both  given  in  [23].  In  this  section,  we  complete  this  collection  by  adding 
another  operator  for  Y  =  T7(curl,  Q^)  intimately  connected  to  the  other  two  operators.  Its 
utility  will  be  clear  in  a  subsequent  section. 

Since  the  Fortin  operators  for  DPG  methods  are  to  be  defined  on  broken  Sobolev  spaces, 
their  construction  can  be  done  focusing  solely  on  one  element.  We  will  now  assume  that 
the  mesh  ■!?/,  is  a  geometrically  conforming  finite  element  mesh  of  tetrahedral  elements.  Let 
Pp{D)  denote  the  set  of  polynomials  of  degree  at  most  p  on  a  domain  D  and  let  Np(D)  = 
Pp-i(D)3  +  x  x  Pp-i(D)3  denote  the  Nedelec  [27]  space.  For  domains  D  c  Mn,  n  =  2,3,  let 
RP(D)  =  Pp-i(D)n  +  xPp-i(D)  denote  the  Raviart-Thomas  [28]  space.  We  use  77p  to  denote 
the  L 2  orthogonal  projection  onto  Pp(K).  From  now  on,  let  us  use  C  to  denote  a  generic 
constant  independent  of  hx  =  diarn  K.  Its  value  at  different  occurrences  may  differ  and  may 
possibly  depend  on  the  shape  regularity  of  K  and  the  polynomial  degree  p. 

Theorem  5.1.  On  any  tetrahedron  K ,  there  are  operators 

■■  h\k)  -  pp+3(k), 

n™l  :  77(curl,  K)  -  Np+3(K), 

Bp+z  :  T7(div,  TT)  — >  i7p+3(7i), 
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such  that  the  norm  estimates 


(5.1a) 

(5.1b) 

(5.1c) 

hold,  the  diagram 
0 - *  H\K)fKL 

grad 


II  np+3  v\\hRK)  <  C \\v\\Hi(K), 
-^p+3-^lli?(curl,^)  ^  C'll^llfl'(curl,/<)) 

II  Hp+3  <Jfl|.ff(div,.K')  ^  C|kll#(div,K)) 


H (curl,  K)  H(dW,K) 


nP+i 


Ijdiv 
UP+ 3 


L\K ) 

Hp+2 


P+3 

0 - *  Pp+3(X)/M 

commutes,  and  these  identities  hold  for  any  v  e  H1(K),  E  e  H( curl,  AT),  and  r  e  i/(div,  AT): 


iVP+3(^)  iW*0  +p+2(A') 


(5.2a)  (<?,  il^a3d  u  -  u)x  =  0 

(5.2b)  (n-a,II^v-v)  =  0 

(5.2c)  (</,  ^7p+3  E  —  E)k  =  t) 

(5. 2d)  <(n  x  A)  x  n,  n  x  (77^  A  -  A)>  =  0 

(5.2e)  (g,i7^v3r-r)A- =  0 

(5.2f)  <n-(ilpd;v3r-r),/i>  =  0 


V?ePH(/f), 

V  n  ■  (7  e  trfvPp+i(X), 

V  <?  e  Pp(K)3, 

V  (n  x  A)  x  n  e  tr^rl  T  Ap+i(A")3, 
VgePP+1(A')3, 
V^etrJ.adPp+2(A'). 


To  provide  a  constructive  proof  of  Theorem  5.1,  we  will  exhibit  Fortin  operators.  We  will 
use  the  exact  sequence  properties  of  the  finite  element  spaces  appearing  as  codomains  of  the 
operators  in  the  theorem.  We  cannot  use  the  canonical  interpolation  operators  in  these  finite 
element  spaces  because  they  do  not  satisfy  (5.1).  Hence  we  will  restrict  the  codomains  of 
our  operators  to  the  following  subspaces  whose  construction  is  motivated  by  zeroing  out  the 
unbounded  degrees  of  freedom. 


Bpr+:j(K)  =  {v  e  Pp+3(K)  :  v  vanishes  on  all  edges  and  vertices  of  A}, 

| Be  Np+z(K)  :  t  ■  E  =  0  on  all  edges  of  K, 

f  fn-  curl  E  =  0  M(f  e  Pp'+2  [dK)  and 
JdK 

I  ((n  x  E)  x  n)  ■  r  =  0  Mr  e  Ri(dK)  1, 

JdK  ) 


BZW  = 


B%(K)  =  jr  e  Rp+3(K )  :  </>n  •  r  =  0  Mfe  P^+2{dK)  j . 

Here  t  denotes  a  tangent  vector  along  the  underlying  edge,  R\(dK)  =  {r  :  r\f  e  R\(f)  for  all 
faces  /  of  dK},  and  P^’+2{dK)  and  P^+2(dK)  are  defined  as  follows.  To  simplify  notation, 
let  Pp(dK)  =  trf((v  R,p+\(K)  (the  space  of  functions  on  dK  that  are  polynomials  of  degree 
at  most  p  on  each  face  of  dK)  and  let  Pp(dK)  =  tr*ad  Ap(A').  Let  Pp’1(dK)  denote  the 
L2{dK)- orthogonal  complement  of  Pp(dK)  +  Po(dK)  in  Pp(dK)  and  let  Pp(dK)  denote  the 
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L2 ( dK) -orthogonal  complement  of  Pp(dK)  in  Pp(dK).  The  following  result  is  proved  in  [23, 
Lemma  3.2], 

Lemma  5.2.  For  any  v  e  Hl(K),  there  is  a  unique  17 qV  in  Bp^(K)  satisfying 
(5.3a)  (q,  IIq  v  —  v)k  =  0  for  all  q  e  Pp-\(K) , 

(5.3b)  (n  ■  a,  77q  v  —  v}  =  0  for  all  n  ■  a  e  tr^v  Rp+i(K),  and 

(5.3c)  ||  LZq  v\\l2(k)  +  hK ||  grad  ilf  u || L2(K)  ^  C  (\\v\\L2(K]  +  hK\\  gradu||L2(fi:))  . 

We  define  /7”™1  as  a  minor  modification  of  the  analogous  operator  in  [23].  Given  any 
v  6  iL^iv),  first  compute  its  mean  value  on  the  boundary 


l  r 

mdK(v )  =  ——  V, 

|  J8K 


then  split  v  =  vq  +  m^(u)  where  vq  =  v  —  m8K(v)  has  zero  mean  trace,  and  finally  define 

(5.4)  i7^a3d  v  =  ilf  u0  +  mdK{v). 

Lemma  5.3.  v  satisfies  (5.2a)-(5.2b)  and  (5.1a)  for  any  v  e  Hl(K ) 

Proof.  Since  IJ^fv  —  v  =  ilg  vo  —  vq,  equations  (5.3a)  and  (5.3b)  immediately  yield  (5.2a) 
and  (5.2b).  To  prove  the  norm  estimate  (5.1a),  note  that  standard  scaling  arguments  imply 

\K\  9  9 

(5.5)  \\mdK(v)\\2L2^  WV\\2L2(dK)Jfif^  ^  ^(IMI L2(K)  + grad  u|||2(ftr)) , 

(5.6)  ||u  -  mdK{v) \\l2{k)  <  ChK\\  gradu||L2(K). 
for  all  v  in  Hl(K).  Combining  (5.5)  and  (5.3c),  we  get 

II  np+3  v\\l2(k )  <  C  (\\v\\L2{K)  +  hK II  gradu||L2W)  , 
while  combining  (5.6)  and  (5.3c), 

hK ||  grad  (77®^  v)\\L2(K)  =  hK\\  grad77f(u  -  m8K{v)) \\l^k),  by  (5.4) 

C  ( || -  m8K{v)\\L2(K)  +  hx ||  grad(u  -  mgK(v))  \\l2(K))  by  (5.3c) 

ChK\\gr&dv\\L2(K)  by  (5.6). 

These  estimates  together  prove  (5.1a).  □ 

The  next  lemma  is  proved  in  [23,  Lemma  3.3].  It  defines  77)^  exactly  as  in  [23]. 

Lemma  5.4.  Any  a  e  Bp'fi:i(K)  satisfying 

(5.7a)  (q,  cr)K  =  0  V  q  e  Pp+l(K)3, 

(5.7b)  <n  ■  a,  fi)  =  0  V/ie  trjad  Pp+2(K), 

vanishes.  Moreover,  for  any  r  e  77(div,  K),  there  is  a  unique  function  I7p+3T  in  Bp'fi?j(K) 
satisfying  (5.2e)-(5.2f).  It  also  satisfies  (5.1c). 

The  remaining  operator  will  be  defined  after  the  next  result.  It  is  modeled  after  the 
previous  two  lemmas,  but  requires  considerably  more  work. 
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Lemma  5.5.  Any  E  e  B™\{K)  satisfying 

(5.8a)  ((j),  E)k  =  0  \ff>  e  Pp(K)3 

(5.8b)  <//,  n  x  E)  =  0  V/je  tr*url  T  Pp+1(K)3, 

vanishes. 


Proof.  Integrating  by  parts  twice  and  using  (5.8b),  we  have 

(5.9)  ifn-  curl  E  =  (curl  E,  grad  iP)k  =  (n  x  P,  grad^)  =  0  V?/>  £  Pp+2(K). 

J8K 

In  addition,  by  Stokes  theorem  applied  to  one  face  /  of  K ,  we  have 

(5.10)  j  nn-cuilE=  j  nE-t  =  0  VkePo^SK), 

Jf  Jdf 

since  E-t  =  0  on  all  edges  by  the  definition  of  B™\(K).  The  definition  of  B™\{K)  also  gives 

(5.11)  f  (fn-  curl  E  =  0  V0  £  P0;^ (dK) . 

J8K 

Since  Pp+2(dK)  =  Pp^(dK)+P£+2(dK)  +  Po(dK),  equations  (5. 9), (5. 10)  and  (5.11)  together 
imply 


(5.12) 


n  ■  curl  E  =  0 


VifePp+2{dK). 


Since  n  ■  curl  E  e  Pp+2(dK),  we  thus  find  that  n  ■  curl E  =  0  on  dK.  This  implies  that  the 
tangential  component  of  E  on  dK ,  namely  ET  =  (n  x  E)  x  n,  has  vanishing  surface  curl,  so 
it  must  equal  a  surface  gradient,  i.e. ,  ET  =  gradT  v  for  some  v  £  Pp+3(dK).  Moreover,  since 
Et  vanishes  on  all  edges,  v  may  be  chosen  to  be  of  the  form  v  =  bjvp  for  some  vp  £  Pp(dK), 
where  hf  is  the  product  of  all  barycentric  coordinates  of  K  that  do  not  vanish  a.e.  on  /. 

To  use  the  remaining  (as  yet  unused)  condition  in  the  definition  of  B™f\{K),  note  that  the 
tangential  component  of  the  coordinate  vector  x,  namely  xT  is  in  R\{dK).  Combining  this 
with  (5.8b),  we  find  that  for  all  ji  £  tr*rl  T  Pp+i(K)3  and  any  n  £  M, 


(5.13) 


where  the  sums  run  over  all  faces  /  of  dK.  For  any  fj,  e  tr*rl  T  Pp+i(P')3,  the  function  /i  x  n 
is  in  the  Raviart-Thomas  space  on  the  closed  manifold  dK  denoted  by  Rp+1(dK).  (Note  that 
unlike  Ri(dK),  this  space  consist  of  functions  with  the  appropriate  compatibility  conditions 
across  edges  of  dK.)  The  surface  divergence  map 

divT  :  Rc+1(dK)  -►  |  w  £  Pp{dK)  :  f  w  =  0 

l  J8K 
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is  surjective.  Hence  the  term  divT(/r  x  n  +  kxt )  appearing  in  (5.13)  spans  all  of  Pp(dK)  as 
p  and  k  are  varied.  Choosing  p  and  k  so  that  divT(^  x  n  +  kxt )  =  vp,  we  conclude  that  vp 
vanishes  and  hence  ET  =  grad T(bfVp)  =  0,  i.e., 

(5.14)  n  x  E  =  0  on  dK. 


Next,  setting  cf  =  curl?’  in  (5.8a)  and  integrating  by  parts,  we  obtain 


(5.15) 


L 


r  ■  curl  E  =  0  Vr  e  Pp+i ( K)c 


From  (5.12)  and  (5.15),  it  follows  that  r  =  curl E  is  in  Bp™s(K),  and  furthermore,  r  satis¬ 
fies  (5.7a)  and  (5.7b).  Hence,  by  Lemma  5.4,  r  vanishes.  Thus  curin'  =  0  and  consequently 
E  =  gradu  for  some  v  e  Pp+s(K).  Furthermore,  by  (5.14),  we  may  choose  v  =  bxvp- 1  for 
some  Vp-\  e  Pp_i(K),  where  bx  is  the  product  of  all  barycentric  coordinates  of  K.  Then 
(5.8a)  implies 

grad(&A"Wp-i)  ■  4>  =  \  bKvp~\  •  div  (j>  —  0  V $  e  Pp{K)3 . 

JK  JK 

It  now  follows  from  the  surjectivity  of  div  :  Pp(K )3  — »  Hp_i(/i)  that  up_i,  and  in  turn 


E  =  grad(6ftmp-i),  vanishes  on  K. 


□ 


The  next  lemma  defines  the  operator  77))” .  It  will  be  useful  to  observe  now  that  for  any 
EeB;f3(K), 

(5.16)  J  ((n  x  E)  x  n)  ■  r  =  0  Vr  e  Ri(dK)  <*=>  |  ((n  x  E)  x  n)  ■  xT  =  0. 

JdK  JdK 

Indeed,  while  the  forward  implication  is  obvious,  the  converse  follows  from  (5.10).  This 
shows  that  the  condition  that  appears  both  in  the  definition  of  B™\(K)  and  in  (5.16)  above, 
actually  amounts  to  just  one  constraint. 


Lemma  5.6.  Given  any  E  e  H(cuil,K),  there  is  a  unique  77^(7]  E  in  ^“^(/i)  satisfy¬ 
ing  (5.2c)-(5.2d). 

Proof.  We  need  to  estimate  n  =  dim Bp"\(K).  First,  note  that  since  Po(dK)  n  Pp+2{dK)  is 
a  one-dimensional  space  of  constant  functions  on  dK, 

dim  Pp+2(dK)  =  dim Pp+2(dK)  -  dim (P£+2(8K)  +  P0(8K )) 

=  dim  Pp+2(dK)  —  dim  Pp+2(dK)  —  dim  Pq(77V)  +  1 

=  6 p  +  11. 

The  tangential  component  E  ■  t  of  any  E  e  Np+s(K )  is  a  polynomial  of  degree  at  most  p  +  2 
on  each  edge,  so  E  ■  t  represents  p  +  3  constraints  per  edge.  Hence,  counting  the  number  of 
constraints  in  the  definition  of  B(fff\{K), 

n  =  dim  B^\(K)  ^  dim Np+3(I<)  -  6 (p  +  3)  -  dim(P^2(dK))  -  1 
=  dim Np+3(K)  -  6 (p  +  3)  -  (6 p  +  11)  -  1 
where  we  have  used  (5.16).  Thus, 

(5.17)  n  ^  dim Np+z(K)  —  12 p  —  30. 
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Next,  we  count  the  number  of  equations  in  (5.2c)— (5. 2d),  namely 
m  =  dim(tr^rl  T  Pp+i(K)3)  +  dim  Pp(K)3 
=  2  dim  Pp+i(dK)  —  6  (p  +  2)  +  dim  Pp(K)3 
=  dim  Np+3(K)  -  Q(p  +  2)  -  6(p  +  3). 


This  together  with  (5.17)  implies  that  m  =  dim  Np+3(K)  —  12p  —  30  ^  n  Thus,  the  sys¬ 
tem  (5.2c)— (5. 2d),  after  using  a  basis,  is  an  m  x  n  matrix  system  of  the  form  Ax  =  d, 
where  x  £  Mn  is  the  vector  of  coefficients  in  a  basis  expansion  of  /If  “3  E  and  d  is  the 
right-hand  side  vector  made  using  the  given  E.  By  Lemma  5.5,  nul(H)  =  0.  Hence 
m  <  n  =  rank(H)  +  nul(H)  =  rank(H)  ^  min(m,  n)  shows  that  m  =  n.  The  system  de¬ 
termining  ZIf“g  E  is  therefore  a  square  invertible  system.  □ 

Lemma  5.7.  For  all  v  e  iL1(iL)/M  and  E  £  //(curl,  K), 

grad(/7^3d  v)  £  B™l3(K),  curl(Z2^  E)  £  B%(K). 

Proof.  Let  e  =  grad(/Iff_ad  v)  =  grad  (/If  vq)  where  /If  vq  £  B^3(K)  is  as  in  (5.4).  Since 
n q  vq  is  constant  along  edges  of  K,  e  must  satisfy  Z  •  e  =  0  along  the  edges.  Moreover, 

(5.18)  m8K{n^v0)  =  0 


due  to  (5.3b).  Hence,  integrating  by  parts  on  on  any  face  /  of  dK , 


l 


xT  =  xT  •  gradT(/lf  u0)  =  -  /If  v0  divT(xT) 
If  If 

Summing  over  all  faces  /  of  dK  and  using  (5.18),  we  conclude  that 


=  -2  /If  ,o. 


Jdj 


eT  ■  x  =  0. 


1PK 
?curl  ( 


Therefore,  to  finish  proving  that  e  e  it  only  remains  to  show  that  K  (fin- curie  =  0 

for  all  <f>  £  Pp’+2(dK).  But  this  is  obvious  from  the  fact  that  e  is  a  gradient. 

Next,  we  need  to  show  that  a  =  curl  (/Iff)]  E)  is  in  Bp™3(K).  Since  it  is  obvious  that 
a  £  Rp+3(K),  it  suffices  to  prove  that 


(5.19) 


f  n  ■  curl  (/If  E)  =  0 

JoK 


for  all  <f>  in  P^~+2(dK).  Note  that  Pp+2(dK)  can  be  orthogonally  decomposed  into  its  subspace 
Po(dK)  n  Pp+2(dK)  and  its  I2(dII)-orthogonal  complement.  The  latter  is  a  subspace  of 
Pp'+2{dK)  where  (5.19)  holds  (since  /If  “3  E  e  B^p\(K)) .  Hence  it  only  remains  to  prove 
that  (5.19)  holds  for  <f  in  Po{dK )  n  Pp+2{dK).  But  Stokes  theorem  shows  that  (5.19)  actually 
holds  for  all  </>  £  P0{dK)  -  cf.  (5.10).  □ 


Lemma  5.8.  For  all  v  £  Hl{K),E  e  //(curl,  K),  and  a  £  //(div,  K), 
(5.20a)  grad  /7f+3d  v  =  ZTf+g  grad  v, 

(5.20b)  curl  E  =  /Id;v3  curl  E, 

(5.20c)  div  /Id+3  cx  =  IIP+ 2  div  a. 
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Proof.  By  Lemma  5.7,  <5i  =  grad(17®+ad  v)  —  grad  v  is  in  We  will  now  show 

that  di  satisfies  (5.8).  Let  4>  e  Pp(K)3  and  consider 

(<t>,  <5i)k  =  {4>,  grad(77^a3d  v  -  v))K  ~  ((/>,  n™ ^(gradu)  -  gradt;)^. 

By  Lemma  5.6,  grad  v  satisfies  (5.2c)-(5.2d),  so  the  last  term  above  vanishes.  Inte¬ 
grating  the  remaining  term  on  the  right-hand  side  by  parts,  and  using  (5.2a)-(5.2b),  we  find 
that 

(5.21)  (0,  grad(17ps™d  v  -  v))K  =  0  V«/»  e  Pp(K)3. 

This  proves  that  (</>,  8\)k  =  0,  i.e.,  (5.8a)  holds.  Next,  for  any  /j  =  tr^rl  T (E) >  E  e  -PP+i(FQ3, 
we  have 

(fi,  n  x  Si)  =  (F  x  n,  grad(17®+ad  v  —  v))  —  (//,  n  x  ( il);"r3 (gr ad  r>)  -  grad  v)). 

The  last  term  vanishes  due  to  (5. 2d).  Moreover, 

(. F  x  n,  grad(ilp)(ad  v  —  v))  =  —(curl  F,  grad(17®+ad  v  —  v))k  =  0 

due  to  (5.21),  so  we  have  proven  that  (5.8b)  holds  as  well.  Hence  by  Lemma  5.5,  S\  =  0. 
This  proves  (5.20a). 

To  prove  (5.20b),  we  proceed  similarly  and  show  that  5 2  =  curl  JT^^E  —  ^p+3  curl  E  is 
zero.  By  Lemma  5.7,  we  know  that  82  £  (K).  so  if  we  prove  that 

(5.22a)  (82,  4>)k  =  0  V(/>  e  Pp+i(^)3, 

(5.22b)  (n  ■  S2,  h)sk  =  0  V/x  £  tr^rad  Pp+2(K). 

then  Lemma  5.4  would  yield  82  =  0.  To  prove  (5.22b), 

<n  •  52,/f)dK  =  (n  •  (curl (i7p+3  E  —  E)  +  (I  —  17d+3)  curl E),/j)8k 

=  (n  •  curl  (IT™1:]  E  -  E),  n)eK  by  (5.7b) 

=  <(77p+3  E  -  E)  x  n,  gradT  fi)8K  =  0  by  (5. 2d). 

To  prove  (5.22a), 

(fc,  4>)k  =  (curl (II™*  E  —  E)  +  (I  —  17^)  curl  E,  <j>)K 

=  (curl (II™1  E-E),  <j>)K  by  (5.2e) 

=  (II™*  E-E,  curl  f)K  +  ((II™*  E-E)xn,<t>)  =  0  by  (5.2c)-(5.2d). 

This  finishes  the  proof  of  (5.22)  and  hence  (5.20b)  follows. 

Finally,  to  prove  (5.20c),  let  £3  =  div  a  —  IIp+2  diver  in  Pp+  2(E).  For  any  w  e 
Pp+2(K),  integrating  by  parts, 

(£3,  w)K  =  (div(77d+3  cr  -  a),w)K 

=  — (77d+v3  a  -  a,  grad  w)K  +  < n  •  (iId+3  a  -  a) ,  w) 

which  vanishes  by  (5.2e)-(5.2f).  Hence  83  =  0  and  (5.20c)  is  proved.  □ 
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Proof  of  Theorem  5.1.  The  lemmas  of  this  section  prove  all  statements  of  Theorem  5.1  ex¬ 
cept  (5.1b).  To  prove  (5.1b),  we  use  a  scaling  argument  and  the  commutativity  properties  of 
Lemma  5.8.  Let  K  denote  the  unit  tetrahedron  and  let  the  IL(curl) -Fortin  operator  on  K, 
defined  as  above,  be  denoted  by  LZ^g  .  By  the  unisolvency  result  of  Lemma  5.6  and  finite 
dimensionality,  there  is  a  C*o  >  0  such  that 

I  rrcurl  771 II 2  ft  II  fp\ |2 

II  11  p+2,  Wl //(curl, if)  '-y°  II  ^  II  //(curl, if) 

for  all  E  G  iL(curl,  K).  Let  Sk  ■  K  — >  K  be  the  one-to-one  affine  map  that  maps  I\  onto 
a  general  tetrahedron  K.  For  E  :  K  — *•  M3,  define  <1>(E)  =  ( S'K)t{E  o  Sk),  and  \\E\\2K  curl  = 
hK2\\E\\2L2iK)  +  ||  curl E 1 1 2l2 ^Ky  The  elementary  proofs  of  the  following  assertions  (i)-(iii)  are 
left  to  the  reader. 

(i)  E  g  B™\(K)  if  and  only  if  §(E)  e  B™"\(k). 

(ii)  There  are  constants  C\ ,  C2  depending  only  on  the  shape  regularity  of  K  (but  not  on 
hK)  such  that  Ci\\E\\2KtCwl  ^  hK\\^(E)\\2kcm]  ^  C2\\E\2KcmV 

(hi)  *(n^E)  =  n^$(E). 

These  three  statements  imply  that 

^11^+3  E\\% curl  <  hK\\n^HE)\\lcurl  ^  hKC0\\*(E)\\%iCail  ^  C0C2\\E\\2Kcml. 

While  this  immediately  gives  the  needed  estimate  for  the  L2-part,  namely 


(5.23) 


n^E\\L,{K) 


^  C\\E\ 


H(cui\,K)  i 


we  need  to  improve  the  estimate  on  the  curl  to  finish  the  proof:  For  this,  we  use  the  commu¬ 
tativity  property 


(5.24) 


curl  77^3  E\\L2(K)  =  ||  ilp+g  curl  E\\L2(K)  by  Lemma  5.8, 
^  C ||  curl  Fill by  Lemma  5.4. 


The  required  estimate  (5.1b)  follows  from  (5.23)  and  (5.24). 


□ 


Before  concluding  this  section,  let  us  illustrate  how  to  use  Theorem  5.1  for  error  analysis 
of  DPG  methods  by  an  example. 

Example  5.9  (Primal  DPG  method  for  the  Dirichlet  problem) .  Consider  the  broken  variational 
problem  of  Example  3.4:  Find  ( u,an )  g  X  =  Hl{fl)  x  H~P2{dQh)  such  that  (4.1)  holds  with 

b{(u,an),y)  =  (grad  u,  grad  y)h  +  (an,y)h,  Y  =  Hl(C2h). 

We  want  to  analyze  the  DPG  method  given  by  (4.3)  with 

Xh  =  {( wh ,  n  ■  fh)  g  X  :  wh\K  G  Pp+1(I<)  and  n  ■  fh\dK  e  trfiv  Rp+i{K) 

for  all  (tetrahedral  mesh  elements)  K  G  1?^}, 

Yh  =  {yh  g  H\Qh)  :  yh\K  e  Pp+:i(K)  for  all  K  e  Qh). 
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We  have  already  shown  in  Example  3.4  that  the  inf-sup  condition  required  for  application  of 
Theorem  4.1  holds  (and  Z  =  {0}).  Hence  to  obtain  optimal  error  estimates  from  Theorem  4.1, 
it  suffices  to  verify  Assumption  3.  We  claim  that  Assumption  3  holds  with  77  =  T/)5)™)1 .  Indeed, 

b{  (■ wh ,  n  ■  Th),y  —  y)  =  (grad  wh,  grad (y  -  77^a3d  y)  )h  +  (n  ■  fh,y  -  77 ^ad  y)h, 

=  -(A wh,  y  -  77^a3d  y)h  +  <n  •  rh  -  y  -  77^ad  y)h  =  0 

by  applying  (5.2a)-(5.2b)  element  by  element.  Note  that  here  we  have  used  the  fact  that 
the  discrete  spaces  have  been  set  so  that  —  Ewh\x  £  Pp-i{K )  and  n  ■  f)x  —  n  ■  grad  Wh  is  a 
polynomial  of  degree  at  most  p  on  each  face  of  dK  (i.e.,  it  is  in  tr^iv  Rp+i(K)),  allowing  us  to 
apply  (5.2a)— (5.2b).  Applying  Theorem  4.1,  we  recover  the  error  estimates  for  this  method, 
originally  proved  in  [17].  /// 


6.  Maxwell  equations 

In  this  section,  we  combine  the  various  tools  developed  in  the  previous  sections  to  analyze 
the  DPG  method  for  a  model  problem  in  time-harmonic  electromagnetic  wave  propagation. 

6.1.  The  cavity  problem.  Consider  a  cavity  17,  an  open  bounded  connected  and  con¬ 
tractible  domain  in  R3,  shielded  from  its  complement  by  a  perfect  electric  conductor  through¬ 
out  its  boundary  dft.  If  all  time  variations  are  harmonic  of  frequency  cj  >  0,  then  Maxwell 
equations  within  the  cavity  reduce  to  these: 

(6.1a)  ~iujpH  +  curl  E  =  0  in  17, 

(6.1b)  —lueE  —  curl  77  =  —  J  in  17, 

(6.1c)  n  x  E  =  0  on  dfl. 

The  functions  E,  77,  J  :  17  — ►  C3  represent  electric  field,  magnetic  field,  and  imposed  current, 
respectively,  and  i  denotes  the  imaginary  unit.  For  simplicity  we  assume  that  the  electro¬ 
magnetic  properties  e  and  p,  are  positive  and  constant  on  each  element  of  the  tetrahedral 
mesh  17/j.  The  number  uj  >  0  denotes  a  fixed  wavenumber.  In  this  section  we  develop  and 
analyze  a  DPG  method  for  (6.1). 

Eliminating  77  from  (6.1a)  and  (6.1b),  we  obtain  the  following  second  order  (non-elliptic) 
equation 

(6.2)  curl  p-1  curl  E  —  u2eE  =  f 

where  /  =  iuiJ .  The  standard  variational  formulation  for  this  problem  is  obtained  by  multi¬ 
plying  (6.2)  by  a  test  function  F  £  77(curl,  17),  integrating  by  parts  and  using  the  boundary 
condition  (6.1c):  Find  E  e  77(curl,  17)  satisfying 

(6.3)  (p-1  curl  E,  curl  F)n  -u2(eE,F)n  =  < f,F ) 

for  any  given  /  £  77(curl,  17)7  It  is  well-known  [25]  that  (6.3)  has  a  unique  solution  for 
every  /  £  77(curl,  17/  whenever  oj  is  not  in  the  countably  infinite  set  E  of  resonances  of  the 
cavity  17.  Throughout  this  section,  we  assume  u  $  E.  This  wellposedness  result  provides  an 
accompanying  stability  estimate,  namely  there  is  a  constant  Cu  >  0  such  that 

(6-4)  l|77||//(cm.i!^)  sS  ^ll/ll_ff(CUri,i?)' 
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for  any  /  £  //(curl,  ft)'  and  E  £  //(curl,  12)  satisfying  (6.3).  Note  that  the  stability  constant 
Cw  may  blow  up  as  u  approaches  a  resonance.  We  continue  to  use  C  to  denote  a  generic 
mesh-independent  constant,  which  in  this  section  may  depend  on  and  e  as  well. 


6.2.  Primal  DPG  method  for  the  cavity  problem.  The  primal  DPG  method  for  the 
cavity  problem  is  obtained  by  breaking  (6.3).  Multiply  (6.2)  by  a  (broken)  test  function 
F  e  //(curl,  fth)  and  integrate  by  parts,  element  by  element,  to  get 

(/G1  curlS,  curl F)h  -  (n  x  /r_1curl E,F)h  -  u2(£E,F)h  =  ( f,F)h . 

Now  set  n  x  H  =  ( iu)~1n  x  curl  E  to  be  an  independent  interface  unknown  which  is  to 

be  found  in  Z/_1/2(div,  fth).  This  leads  to  the  variational  problem  (3.4)  with  the  following 
spaces  and  forms: 

(6.5a)  A'o  =  //(curl,  12),  Y  =  //(curl,  fth), 

(6.5b)  X  =  H~1/2( div,  dft h),  To  =  //(curl,  ft), 

(6.5c)  bo(E,F)  =  (/u~1  curl  E,  curl  F)h  -  u2(eE,  F)h, 

(6.5d)  b(H^,F)  =  -iu(H^F)h. 

This  is  the  primal  DPG  formulation  for  the  Maxwell  cavity  problem. 

The  numerical  method  discretizes  the  above  variational  problem  using  subspaces  Xh  cr 
X  =  Xq  x  X  and  Y/,  cr  Y  defined  by 

(6.6a)  Xh  =  {( Eh,n  x  Hh)  e  //(curl,  12)  x  Z/-1/2(div,  dfth)  :  n  x  Hh \dK  £  tr^rl  H  Pp+1(K)3, 

and  E/-, \ k  £  Pp(K)3  for  all  K  e  fth}, 

(6.6b)  Yh  =  {Fh  £  //(curl,  fth)  :  Fh \K  e  Np+3(K)  for  all  K  e  fth}. 

We  have  the  following  error  bound  for  the  numerical  solution  in  terms  of  the  mesh  size 
h  =  max-Ke(}h  hx  and  polynomial  degree  p  ^  1. 

Corollary  6.1.  Suppose  ( Eh,n  x  Hh)  £  Xh  is  the  DPG  solution  given  by  (4.3)  with  forms 
and  spaces  set  by  (6.5)  and  let  ( E,n  x  H)  e  X  be  the  exact  solution  of  (4.1).  Then,  there  is 
a  C  depending  on  oj,  p,  and  the  shape  regularity  of  the  mesh  (but  independent  of  h)  such  that 

\\E  ~  /'ft II //(curl, fi)  +  lln  X  (//  —  Z/h)||tf-l/2(div,dfth) 

^  Chip  [\E\hp+i^)  +  |  cur\E\HP+i(ty  +  |  curl// 1 hp+1{Q) )  • 


Proof.  To  apply  Theorem  4.1,  we  must  verify  the  inf-sup  condition  (4.5)  for  the  broken  form. 
As  in  the  previous  examples,  as  a  first  step,  we  verify  the  inf-sup  condition  for  the  unbroken 
form  stated  in  Assumption  1.  Given  any  E  e  //(curl,  12),  let  fp  e  //(curl,  ft)'  be  defined  by 
(fE,F)  =  (p~1  curl//,  curl F)q  —  u2(eE,  F)q  for  all  F  £  //(curl,  12).  Then,  (6.2)  and  (6.4) 


implies 


Z' lli/(curl,i?)  ^  Co)  1 1 Ie  1 1  fj (curI  fty 


Cu  sup 

FeH{c  url,l7) 


ME,F)\  ^ 

|/1||i/(curl,l7) 


i.e.,  Assumption  1  holds  with  co  =  Cff1.  Assumption  2,  with  c  =  u;-1  is  immediately  verified 
by  (2.11c)  and  (2.10c)  of  Theorem  2.3.  Hence  Theorem  3.1  verifies  (4.5)  and  also  shows  that 
Z  =  {0}.  The  only  remaining  condition  to  verify  before  applying  Theorem  4.1  is  Assump¬ 
tion  3,  which  immediately  follows  by  the  choice  of  spaces  and  Theorem  5.1. 
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Applying  Theorem  4.1,  we  find  that 

II ~  Eh\\2H(curl,Si)  +  \\n  X  (H  -  Hh)YH-i/2(Aivjnh) 

^  C  inl  [ll^  “  ^/illi/(curl,fi)  +  II n  x  H  -  n  X  ^h||^_i/2(di  5J7  ) 
( Gh,nxRh)eXh  L 


Now,  H  =  (iw//)_1  curl/?  is  an  extension  to  i?  of  the  exact  interface  solution  nxH.  Moreover, 
the  interface  function  n  x  Rh  appearing  above  can  be  extended  into  X^h  =  {r  £  //(curl.  1?)  : 
Vhlx  £  F’p+i(K)3J.  Since  the  interface  norm  is  the  minimum  over  all  extensions, 

WE  ~  Eh\\H(cuT\,n)  +  lln  X  (H  -  ^/i)||ff-i/2(diV;5J7a) 


<  C 


inf 

GheXlh 


E-Gh 


2 

H  (cml,f2) 


inf 

RheXSX1 


H  —  Rh 


2 

H  (curl,  f2) 


iC  2 

KeQh 


j  2(si+l) 
nK 


\E\ 


2 


+  curl  E\~Hsl+1(Kj  + 


u  2(s2  +  1) 
nK 


H 


2 

HS2+1(K) 


h2s 2 


curl  H 


2 

Lf/s2+!(A:) 


by  standard  approximation  estimates  (see  e.g.,  [19,  Theorem  8.1])  where  1/2  <  si  ^  p  and 
1/2  <  S2  <  p  +  1.  Hence  the  corollary  follows.  □ 


Remark  6.2.  Unlike  the  standard  finite  element  method,  for  the  DPG  method,  there  is  no 
need  for  h  to  be  “sufficiently  small”  to  assert  a  convergence  estimate  as  in  Corollary  6.1. 


6.3.  Alternative  formulations  of  the  same  problem.  In  Example  3.5,  we  saw  that 
a  single  diffusion-convection-reaction  equation  admits  various  different  formulations.  The 
situation  is  similar  with  Maxwell  equations.  First,  let  us  write  (6.1)  in  operator  form  using 
an  operator  A  (analogous  to  the  one  in  (3.10),  but  now)  defined  by 

iujuH  —  curl  E 
uaeE  +  curl  El 

as  A(H ,  E)  =  (0,  J)  for  some  given  J  in  L2(l?)3.  However,  we  will  not  restrict  to  right-hand 
sides  of  this  form  as  we  will  need  to  allow  the  most  general  data  possible  in  the  ensuing 
wellposedness  studies. 

We  view  A  as  an  unbounded  closed  operator  on  L2(i 2)6  whose  domain  is 
dorri(A)  =  {(//,  E)  £  //(curl,  l?)2  :  n  x  E  =  0  on  dQ}. 


H 

lUJf-L 

—  curl 

H 

E 

curl 

iuje 

E 

It  is  easy  to  show  that  its  adjoint  (in  the  sense  of  closed  operators)  is  the  closed  operator  A* 
given  by 


H 

—lujfr 

curl 

H 

— vjjfiH  +  curl  E 

E 

—  curl 

—iuje 

E 

—lOjeE  —  curl  H 

whose  domain  is  the  following  subspace  of  L2(i?)6: 

dom(A*)  =  {(//,  E)  £  //(curl,  17)2  :  Er  =  0  on  dl2}. 


Classical  arguments  show  that  both  A  and  A*  are  injective.  To  facilitate  comparison,  we  list 
all  our  formulations  at  once,  including  the  already  studied  primal  form. 


27 


Strong  form:  Let  x  =  ( H ,  E )  be  a  group  variable.  Set 

W0  =  H (curl,  S2)  x  H{ curl,  i7),  Y  =  Y0  =  L2(l7)6. 

Note  that  Xq  =  {dom(A),||  •  ||#(curi,j?)}>  i-e.,  dom(A)  considered  as  a  subspace  of 
//( curl,  l?)2  (rather  than  as  a  subspace  of  L2(l?)6).  The  Maxwell  problem  is  to  find 
x  6  Xq.  given  /  e  Yq,  such  that  Ax  =  f.  This  fits  into  our  variational  framework  (3.2) 
by  setting  bo  to 

bo(x,y)  =  (Ax,y)n- 

Primal  form  for  E:  This  is  the  same  as  in  (6.5),  i.e.,  with  the  spaces  as  set  there, 
with  b  set  to  bE{H H,  F)  =  —iuj(H^  ,  F}h  and  bo  set  to 

bE (E,  F)  =  (|U_1  curl  E,  curl  F)q  -  uj2(eE,F)a, 

the  electric  primal  formulation  is  (3.2)  and  its  broken  version  is  (3.4). 

Primal  form  for  H:  Eliminating  E  from  (6.1),  we  obtain  curl  Y  1  curl  II  —  u'2 (iH  = 
curie-1,/  and  a  (possibly  nonhomogeneous)  boundary  condition  on  n  x  e_1  curl H. 
With  this  in  place  of  (6.3)  as  the  starting  point  and  repeating  the  derivation  that  led 
to  (6.5),  we  obtain  the  following  magnetic  primal  form.  Set 

X0  =  H (curl,  f2),  Yo  =  Xo, 

X  =  i/-1/2(div,  f2h),  Y  =  H (curl,  Qh), 

bo  (H,  F)  =  (e-1  curl H,  curl F)h  —  F)h,  bH{E^F)  =  wj(EH,F)h. 

With  b  set  to  bH  and  bo  set  to  bE ,  the  magnetic  primal  formulation  is  (3.2)  and  its 
broken  version  is  (3.4). 

Ultraweak  form:  This  form  is  obtained  by  integrating  by  parts  all  equations  of  the 
strong  form.  Using  group  variables  x  =  ( H ,  E),  y  =  (R,  S )  and  x  =  (HT ,  Er),  set 

(6.7a)  W0  =  L2(l?)6,  Y0  =  //(curl,  Q)  x  //(curl,  f2), 

(6.7b)  Y  =  //(curl,  l?/,)2,  X  =  //_1/2(curl,  dQfi)  x  //_1^2(curl,  dfih), 

(6.7c)  b% (x,  y )  =  (x,  A*y)h,  bu(x,  y)  =  < HT,n  x  S)h  -  (ET,  n  x  R)h 

and  consider  formulations  (3.2)  and  (3.4)  with  bo  =  /[/  and  b  =  bu .  Note  that  in 
the  definition  of  bE,  the  operator  A *  is  applied  element  by  element,  per  our  tacit 
conventions  when  using  the  (•,  ^-notation. 

Dual  Mixed  form:  Among  the  two  equations  in  the  strong  form,  if  one  weakly  imposes 
(by  integrating  by  parts)  the  first  equation  and  strongly  imposes  the  second,  then  we 
get  the  following  dual  mixed  form.  Set 

W0  =  //(curl,  S2)  x  L2(f2)3,  Y0  =  W0, 

W  =  H~1/2( curl,  dnh),  Y  =  //(curl,  S2h)  x  L2(i 7)3 

and  consider  (3.2)  with  bo  =  bE, 

bo  ( (H,  E),  (R,  S) )  =  (udiiH,  R)q  —  (E,mAR)h 
+  (lixeE  +  curl  ZZ,  S)h- 
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Its  broken  version  is  (3.4)  with  b  set  to  bD(ET,  (R,  s) )  =  (ET,n  x  R)^. 

Mixed  form:  Reversing  the  roles  above  and  weakly  imposing  the  second  equation  while 
strongly  imposing  the  first,  we  get  another  mixed  formulation.  Set 

X0  =  L2(C)3  x  H (curl,  12),  Y0  =  X0, 

X  =  H-V2( curl,  dQh),  Y  =  L2(C)3  x  //(curl,  S2h) 

and  consider  (3.2)  with  bo  =  , 

b^(  (//,  E),  (R,  S) )  =  (iupH  -  curl  E,  R)h 

+  ( iojeE,S)n  +  (H,cui\S)h- 

Its  broken  version  is  (3.4)  with  b  set  to  bM(HT ,  (/?,  S ) )  =  (i/T,  n  x  S'/^. 

These  form  a  total  of  six  unbroken  and  five  broken  formulations,  counting  the  already  dis¬ 
cussed  broken  and  unbroken  electric  primal  formulation.  To  analyze  the  remaining  formu¬ 
lations,  let  us  begin  by  verifying  Assumption  1  for  all  the  unbroken  formulations.  To  this 
end,  label  the  statement  of  Assumption  1  with  bo  set  to  the  above-defined  6q  as  “(2)”  for  all 
I  e  {E,  H,  S,U,  D,  Mj.  Then  (cf.  Figure  1)  we  have  equivalence  of  statements  (D),  (E), . . . 
as  proved  next. 

Theorem  6.3.  The  following  implications  hold: 


Proof.  We  begin  with  the  most  substantial  of  all  the  implications,  which  allows  us  to  go  from 
the  strongest  to  the  weakest  formulation.  When  there  can  be  no  confusion,  let  us  abbreviate 
Cartesian  products  of  L2(i?)  as  simply  L  and  write  W  for  //(curl,  i?)  x  //(curl,  f2).  Clearly,  W 
is  complete  in  the  //(curl,  f?)2-norm.  It  is  easy  to  see  that  the  graph  norms  (||x||2  +  ||At||2  )1//2 
and  (|i.xj|2  +  ||A*x|[2  )x/2  are  both  equivalent  to  the  //(curl,  f?)2-norm,  so  W  is  a  Hilbert  space 
in  any  of  these  norms.  These  norm  equivalences  show  that  the  inf-sup  condition  ( S )  holds  if 
and  only  if 

(6.8)  CIMIl  <  1 1 -Ax 1 1 l  V.x  6  dorn(A). 

(5)  =>  (U):  The  bound  (6.8)  implied  by  (5)  shows  that  the  range  of  A  is  closed.  By  the 
closed  range  theorem  for  closed  operators,  range  of  A*  is  closed.  Since  A*  is  also  injective,  it 
follows  that 


Vy  e  IT. 
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Thus,  to  complete  the  proof  of  (U),  it  suffices  to  show  that 


(6.10) 


inf  sup 

xeL  yeW 


I  (x,A*y)n 

IMUIMIw 


inf  sup 

xeL 


\(x,A*y)n\ 
IMIlIMIw  ' 


For  completeness,  we  now  describe  the  standard  argument  that  shows  that  one  may  reverse 
the  order  of  inf  and  sup  to  prove  (6.10).  Viewing  A*  :  W  — >  L  as  a  bounded  linear  operator,  we 
know  that  it  is  a  bijection  because  of  (6.8)  and  (6.9).  Hence  (H*)-1  :  L  — >  W  is  bounded.  The 
right-hand  side  of  (6.10)  equals  its  operator  norm  ||(H*)^1||.  The  left  hand  side  of  (6.10)  equals 
the  operator  norm  of  the  dual  of  (H*)-1  (considered  as  the  dual  operator  of  a  continuous 
linear  operator  with  L  as  the  pivot  space  identified  to  be  the  same  as  its  dual  space).  The 
norms  of  a  continuous  linear  operator  and  its  dual  are  equal,  so  (6.10)  follows. 

(U)  =>  ( D ):  Let  x  =  ( H ,  E)  and  y  =  ( R ,  S)  be  in  W  =  //(curl,  1?)  x  //(curl,  17).  Clearly, 
W  is  contained  in  Y{f  =  //(curl,  17)  x  L2(l 7)3.  Because  of  the  extra  regularity  of  S,  we  may 
integrate  by  parts  the  last  term  in  the  definition  of  6^((//,  E),  ( R ,  5))  to  get  that  b q  (x,  y)  = 
b(l(x,y).  Hence  using  (U), 


(6.11) 


sup 

y£Y0D 


K(x,y)\ 

\\y\\yn 


> 


sup 

yeW 


K(x,y)\ 

h\\w 


sup 

yeW 


\bp(x,y)\ 

\\y\\w 


>C\\x\\L. 


Thus,  to  finish  the  proof  of  (Z7),  we  only  need  to  control  curl//  using  the  last  term  of  b^. 

II  r„rl  mi  KUH.E),  (0 ,S))  -  (mieE, S) 

SeL  11*8 \\L  SeL  11*8  \\l 

^  sup 
ycY0D 

Using  (6.11)  to  bound  the  last  term,  the  proof  of  ( D )  is  finished. 

(. D )  ==>  (//):  For  any  H  e  //(curl,  17),  set  Ih(R)  =  (e_1  curl//,  curl/?) q  —  cu2(/r/Z,  R)q. 
We  need  to  prove  (//),  which  is  equivalent  to 


\bp(x,y)\ 

112/11^0° 


WiujeE  II 


(6-12)  \\H\\H(curi,n)  ^  C\\£H\\H(cuii,ay- 

Introducing  a  new  variable  E  =  — (zcue)-1  curl  //,  we  find  that  ( iuj)~1Ih(R )  =  ~(E,  curl  R)a  + 
R)q.  Hence 

E),  (R,  S ))  =  (*u;)”1/j/(//)  V(/2,  S )  e  //(curl,  17)  x  L2(C)3. 

Hence  (6.12)  immediately  follows  from  ( D ). 

(//)  =>  ( S ):  To  prove  the  inf-sup  condition  (5),  it  is  enough  to  prove  (6.8)  for  all  x  = 
(//,  E)  e  W.  Given  any  F,G  e  L2(l7)3,  the  equation  Ax  =  (F,  G)  is  the  same  as  the  system 


(6.13a)  KjJ/j,H  —  curl  E  =  F, 

(6.13b)  iweE  +  curl  H  =  G. 

We  multiply  (6.13b)  by  the  conjugate  of  e_1curl//,  for  some  R  e  //(curl,  17)  and  integrate 
by  parts,  while  we  multiply  (6.13a)  by  — tuiR  and  solely  integrate.  The  result  is 

(e-1  curl  //,  curl  R)n  +  (u o  curl  E,R)q  =  (G,e  1  curl  R)q, 

-(u2ixH,  R)n  -  (zu  curl E,R)n  =  ~(F,iuR)a- 


30 


C.  CARSTENSEN,  L.  DEMKOWICZ,  AND  J.  GOPALAKRISHNAN 


Adding  the  above  two  equations  together,  we  get  the  primal  form  6j^(#,  #)  =  £(R)  where 
£(R)  =  (G,  £_1  curl  i?)  —  (F,tuiR)s ?•  Hence  the  given  inf-sup  condition  (#)  implies 


C\\H\\Hicur\,n)  <  sup 

ReH  (curl,/?) 


K(h,r)\ 


||#| 


H(  curl,!?) 


sup 

i?ei?(curl,G) 


\m\ 

l|-R||.ff(curlIJ?) 


Since  \£(R)\  ^  CdjFlIx,  +  ||Gf||L)||-R||_ff(curi,j7),  this  provides  the  required  bound  for  ||#||£. 
Since  curl#  is  also  bounded,  equation  (6.13b)  yields  a  bound  for  \\E\\l-  Combining  these 
bounds,  (6.8)  follows. 

To  conclude  the  proof  of  the  theorem,  we  note  that  the  proofs  of  the  implications  ( U )  => 
(M),  (M)  =>  (S),  (E)  =>  ( S )  are  similar  to  the  proofs  of  ( U )  =>  (D),  (D)  =>  (H), 
and  (#)  =>  (S'),  respectively.  □ 


Theorem  6.3  verifies  Assumption  1  for  all  the  formulations  because  we  know  from  (6.4) 
that  (E)  holds.  Assumption  2  can  be  easily  verified  for  all  the  broken  formulations  using 
Theorem  2.3.  Assumption  3  can  be  verified  using  Theorem  5.1.  Hence  convergence  rate 
estimates  like  in  Corollary  6.1  can  be  derived  for  each  of  the  broken  formulations.  We  omit 
the  repetitive  details. 


7.  Numerical  studies 

In  this  section,  we  present  some  numerical  studies  focusing  on  the  Maxwell  example.  Nu¬ 
merical  results  for  other  examples,  including  the  diffusion-convection-reaction  example,  can 
be  found  elsewhere  [10,  13].  The  numerical  studies  are  not  aimed  at  verifying  the  already 
proved  convergence  results,  but  rather  at  investigations  of  the  performance  of  the  DPG 
method  beyond  the  limited  range  of  applicability  permitted  by  the  theorems.  All  numerical 
examples  presented  in  this  section  have  been  obtained  with  hp3d,  a  3D  finite  element  code 
supporting  anisotropic  h  and  p  refinements  and  solution  of  multiphysics  problems  involv¬ 
ing  variables  discretized  compatibly  with  the  H 1  (17)-# (curl,  l?)-#(div,  1?)  exact  sequence  of 
spaces.  The  code  has  recently  been  equipped  with  a  complete  family  of  orientation  embedded 
shape  functions  for  elements  of  many  shapes  [21].  The  remainder  of  this  section  is  divided 
into  results  from  two  numerical  examples. 

Example  7.1  (Smooth  solution).  We  numerically  solve  the  time- harmonic  Maxwell  equations 
setting  material  data  to 

e  =  n  =  1,  oj  =  1, 

and  1?  to  the  unit  cube.  To  obtain  the  unit  cube  was  partitioned  first  into  five  tetrahedra: 
four  similar  ones  adjacent  to  the  faces  of  the  cube,  and  a  fifth  inside  of  the  cube.  We  have 
used  the  refinement  strategy  of  [24]  to  generate  a  sequence  of  successive  uniform  refinements. 
On  these  meshes,  consider  the  primal  DPG  method  for  E,  described  by  (6.5),  with  data  set 
so  that  the  exact  solution  is  the  following  smooth  function. 


Ei  =  sin7rxi  sin7TX2  sin  71x3,  E2  =  #3  =  0  . 
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(a)  Tetrahedral  meshes. 


Number  of  dof  (N) 


(b)  Hexahedral  meshes. 

Figure  2.  Rates  for  the  case  of  smooth  solution  and  primal  formulation. 


Instead  of  the  pair  of  discrete  spaces  (6.6)  that  we  know  is  guaranteed  to  work  by  our 
theoretical  results,  we  experiment  with  these  discrete  spaces: 

(7.1a)  Xh  =  {(Eh,  n  x  Hh)  e  (curl,  12)  x  tf-1/2(div,  df2h)  :  n  x  Hh \8K  e  tr^lH  Np(K), 

and  Eh\x  £  Np(K)  for  all  K  e  17 /J, 

(7.1b)  Yh  =  {Fh  e  i7(curl,  Qh)  :  Fh \K  e  1Vp+2(/F)  for  all  77  e  Qh). 
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(a)  Tetrahedral  meshes. 


Number  of  dof  (N) 


(b)  Hexahedral  meshes. 

Figure  3.  Rates  for  the  case  of  smooth  solution  and  ultraweak  formulation. 


The  observed  rates  of  convergence  of  the  error  \\E  —  Eh\\H(cmi,oh)  and  the  residual  norm  r/ 
are  shown  in  Figure  3a.  The  rates  are  optimal.  This  suggests  that  the  results  of  Corollary  6.1 
may  hold  with  other  choices  of  spaces.  The  problem  of  finding  the  minimal  DPG  test  space  for 
which  optimal  convergence  rates  can  be  obtained  for  the  Maxwell  problem  remains  unsolved. 

We  also  present  similar  results  obtained  using  cubic  meshes  using  curl,  h?)-conforming 
Nedelec  hexahedron  of  the  first  type.  Namely,  Xh  and  Y h  are  set  by  (7.1)  after  revising  Np(K) 
to  x  QPtP_i)P(K)  x  QPtP:P_i(K)  where  Qitm,n(K)  denotes  the  set  of  polynomials 
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Figure  4.  Construction  of  Fichera  oven. 

of  degree  at  most  l,  rn,  and  n  in  the  xi,x2  and  X3  directions,  respectively.  The  convergence 
rates  reported  in  Figure  3b  are  again  optimal. 

Before  concluding  this  example,  we  also  report  convergence  rates  obtained  from  the  ultra- 
weak  formulation  of  (6.7).  The  discrete  spaces  are  now  set  by 

(7.2a)  Xh  =  {(E,H,ET,HT)  e  L2(C)3  x  L2(f2)3  x  i7~1/2(curl,  d!2h)  x  #~1/2(curl,  dS2h)  : 

E\k,  H\k  e  PP-i(Kf,  Et \dK,  HT\dK  e  tr^rl>T  Np(K)  for  all  K  e  Qh), 
(7.2b)  Yh  =  {(F,  G)  e  H{ curl,  Qh)  x  H{ curl,  S2h)  :  F\K,  G\K  e  Np+2(K )  for  all  K  e  Qh}. 

Recall  that  the  DPG  computations  require  a  specification  of  the  Y -norm.  Using  the  obser¬ 
vation  (made  in  the  proof  of  Theorem  6.3)  that  the  adjoint  graph  norm  is  equivalent  to  the 
natural  norm  in  H( curl,  C)2,  we  set 

\\(E,H)\\y  =  ^  ( \\e\\2l2({2 )  +  II^IIl2(U)  +  -  curlUjl^^)  +  \\tueE  +  curli7|||2(fi)^) 

KeQh 

in  all  computations  involving  the  ultraweak  formulation.  The  results  reported  in  Figure  3 
again  show  optimal  convergence  rates.  Note  that  only  the  errors  in  the  interior  variables 
E  and  H  (in  L2(h?)-norm)  are  reported  in  the  figure.  To  compute  errors  in  the  interface 
variables,  we  must  compute  approximations  to  fractional  norms  carefully  (see  [7]  for  such 
computations  in  two  dimensions).  Since  the  code  does  not  yet  have  this  capability  in  three 
dimensions,  we  have  not  reported  the  errors  in  interface  variables.  /// 

Example  7.2  (Singular  solution).  To  illustrate  adaptive  possibilities  of  DPG  method  and  the 
difference  between  different  variational  formulations,  we  now  present  results  from  a  “Fichera 
oven”  problem.  We  start  with  with  the  standard  domain  with  a  Fichera  corner  obtained 
by  refining  a  cube  (0,2)3  into  eight  congruent  cubes  and  removing  one  of  them.  We  then 
attach  an  infinite  waveguide  to  the  top  of  the  oven  and  truncate  it  at  a  unit  distance  from 
the  Fichera  corner,  as  shown  in  Figure  4.  Setting 


e  =  n  =  1,  u>  =  5 
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Figure  5.  Every  other  iterate  (iterates  1,3, 5, 7, 9  and  11)  from  adaptive  algo¬ 
rithm  applied  to  solve  the  Fichera  oven  problem  with  the  primal  formulation. 
Meshes  (left)  and  the  corresponding  real  part  of  E\  are  shown. 
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(a)  Adaptivity  for  primal  formulation.  (b)  Adaptivity  for  ultraweak  formulation. 

Figure  6.  Convergence  of  residual  during  adaptive  iterations  for  the  Fichera  oven. 

we  drive  the  problem  with  the  first  propagating  waveguide  mode, 

Ei  =  sin  7rx'2,  E 2  =  E%  =  0 

which  is  used  for  non-homogenous  electric  boundary  condition  n  x  E  =  n  x  ED  across  the 
waveguide  section.  Analogous  to  a  microwave  oven  model,  we  set  the  homogeneous  perfect 
electric  boundary  condtion  n  x  E  =  0  everywhere  else  on  the  boundary.  The  above  material 
data  correspond  to  about  0.8  wavelengths  per  unit  domain.  In  all  the  reported  computations, 
we  start  with  a  uniform  mesh  of  eight  quadratic  elements  that  clearly  does  not  even  meet 
the  Nyquist  criterion.  We  expect  the  solution  to  develop  strong  singularities  at  the  reentrant 
corner  and  edges,  but  we  do  not  know  the  exact  solution. 

First,  we  report  the  results  from  the  electric  primal  formulation,  choosing  spaces  again 
as  in  (7.1).  Figure  5  presents  the  evolution  of  the  mesh  along  with  the  corresponding  real 
part  of  the  first  component  of  electric  field  E\.  Since  we  do  not  have  the  exact  solution  for 
this  problem,  we  display  convergence  history  using  a  plot  of  the  evolution  of  the  computed 
residual  rj  in  Figure  6a.  (Recall  that  theoretical  guidance  on  the  similarity  of  behaviors  of 
error  estimator  rj  and  the  error  is  provided  by  Theorem  4.1.)  Clearly,  the  figure  shows  the 
residual  is  being  driven  to  zero  during  the  adaptive  iteration. 

Next,  we  solve  the  same  problem  using  the  ultraweak  formulation  with  the  spaces  set  as 
in  (7.2)  and  the  Y -norm  set  to  the  adjoint  graph  norm  as  in  the  previous  example.  The 
convergence  history  of  the  residual  norm  ?y  is  displayed  in  Figure  6b.  The  evolution  of  the 
mesh  along  with  the  real  part  of  E\  is  illustrated  in  Figure  7. 

It  is  illustrative  to  visualize  the  difference  between  the  two  different  DPG  formulations  and 
the  accompanying  convergence  in  different  norms.  Figure  8  presents  a  side-by-side  comparison 
of  the  real  part  of  the  electric  field  component  E\  obtained  using  the  primal  (left)  and 
ultraweak  (right)  formulations.  The  same  color  scale  (min  =  —  1,  max  =  1)  is  applied  to  both 
solutions  in  this  figure  (whereas  the  scales  of  Figures  5  and  7  are  not  identical).  Obviously,  the 
meshes  are  different,  but  they  are  of  comparable  size,  so  we  believe  the  comparison  is  fair.  The 
primal  method,  which  delivers  solution  converging  in  the  stronger  curl,  i7)-norm,  “grows” 
the  unknown  solution  slower,  whereas  the  ultraweak  formulation  converging  in  the  weaker 
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Figure  7.  Iterates  1,3, 5, 7,  and  9  from  adaptive  algorithm  for  Fichera  oven 
problem  with  the  ultraweak  formulation.  Meshes  (left)  and  the  corresponding 
real  part  of  E\  are  shown. 


L2(I?)-norm  seems  to  capture  the  same  solution  features  faster.  Both  methods  ultimately 
deliver  the  same  solution  but  at  different  speeds  and  the  ultraweak  formulation  seems  to  be 
a  winner.  Recall  that  the  number  of  interface  unknowns  for  both  formulations  is  identical, 
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(a)  Primal  (left)  and  ultraweak  (right)  iterate  6. 


(c)  Primal  (left)  and  ultraweak  (right)  iterate  8.  (d)  Primal  (left)  and  ultraweak  (right)  iterate  9 

Figure  8.  Same  scale  comparison  of  real  part  of  E\  computed  by  the  primal 
and  the  ultraweak  formulations.  A  middle  slice  (through  singular  vertex)  of 
adaptive  iterates  6,7,8  and  9  are  shown. 


but  the  total  number  of  unknowns  for  the  ultraweak  formulation  is  higher,  i.e. ,  the  ultraweak 
formulation  requires  a  larger  number  of  local  (element-by-element)  computations.  /// 
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